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Research  on  Deterministic  Methods  of  Seismic  Source  Identification 

Report  No.  1  -  Summary 


In  order  to  interpret  seismic  event  discrimination  in  terms  of  the 
physical  properties  of  the  source  and  to  be  able  to  establish  new  discrimin¬ 
ation  techniques  we  have  generalized  seismic  source  models  based  on  relaxation 
source  theory  to  include  the  effects  of  non-homogeneous  initial  prestress. 

In  particular  we  have  considered  the  effects  of  strongly  concentrated  prestress 
in  the  vicinity  of  the  shatter  zone  produced  by  an  explosion. 

The  important  result  from  the  work  so  far  completed  are:  (1)  the  spectra 
of  P  and  S  waves  radiated  due  to  stress  relaxation  effects  can  be  strongly 
peaked,  with  the  nature  of  the  peaking  being  azimuthally  dependent  in  general 
and  quite  strongly  dependent  on  the  size  and  location  of  the  initial  stress 
concentration;  (2)  the  comers  or  peak  frequencies  of  the  P  and  S  wave 
radiation  are  different  from  one  another  (S  lower)  and  are  both  shifted  to 
higher  frequencies  when  the  stress  concentration  is  close  to  the  shatter  zone. 

In  addition,  the  corner  or  peak  frequency  value  is  related  to  the  size  of  the 
stress  concentration  rather  than  to  the  size  of  the  shatter  zone;  (3)  the 
pattern  of  first  motions  from  the  tectonic  release,  when  the  prestress  is 
inhomogeneous,  is  not  pure  quadrupole  with  higher  order  multiples  also  involved. 
The  ordinary  quadrupolar  pattern  predicted  from  a  homogeneous  prestress  can 
be  strongly  distorted  when  the  prestress  is  concentrated  and  can  be  highly 
non-quadrupole  in  form. 

Considering  both  the  explosion  generated  direct  compressional  wave 
field  and  the  field  produced  by  such  tectonic  effects  together,  the  consequences 
for  discrimination  and  explosion  yield  determination  are: 

(1)  Short  period  perturbations  in  the  wave  train  can  be  expected  to  be 
very  complex  due  to  dependence  on  stress  concentration  effects.  However, 
the  perturbations  should  be  small  to  moderate  for  the  first  cycle  of 

the  P  wave  motion,  while  being  significantly  larger  for  the  later  part 
of  the  P  wave  train.  Body  wave  magnitude  measured  from  the  first  P 
wave  cycle  should,  therefore,  be  minimally  perturbed  by  stress  relaxation 
effects. 

(2)  Long  period  surface  wave  radiation  can  be  strongly  perturbed  by 
tectonic  release  effects  within  the  whole  measureable  low  frequency 

band  (i.e.g,  from  approximately  5  sec  to  60  sec  in  period).  The  perturba¬ 
tions  in  the  observed  Rayleigh  wave  forms  can  be  such  as  to  add  to, 
or  subtract  from,  the  explosive  generated  Rayleigh  wave  depending  on  the 
orientation  and  magnitude  of  the  prestress  in  the  vicinity  of  the  explosion. 
The  magnitude  of  this  effect  can  be  very  (unacceptably)  large.  In  those 
cases  where  Love  waves  are  significant,  so  that  tectonic  release  is 
involved,  then  yield  estimation  using  M  can  only  be  made  after  correction 
of  the  Rayleigh  wave  measurement  using  the  observed  Love  wave  to  deduce 
the  size  and  configuration  of  the  tectonic  source. 

Predictions  of  the  radiated  seismic  wave  field  from  both  explosion  and 
earthquake  models  in  the  regional  and  teleseismic  distance  ranges,  for 
layered,  anelastic  earth  models,  have  been  used  to  predict  entire  synthetic 
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seismograms  for  purposes  of  comparisons  with  observations.  Both  full  wave 
asymptotic  and  locked  mode  approximation  methods  are  being  used.  Examples  of 
synthetics  and  comparisons  to  observations  in  the  regional  distance  range 
suggest  that  most  of  the  seismogram  can  be  understood  using  current  source 
and  medium  models.  Some  specific  applications  of  this  capability  have  been 
directed  toward  providing  an  understanding  of  the  large  amplitude  Pg  signals 
in  the  regional  distance  range,  and  for  studies  of  regional  anleastic  char¬ 
acteristics  of  the  earth.  The  results  of  these  two  applications  have  been 
to  show  that:  Pg  is  comprised  of  a  large  number  of  high  modes,  which  can  also 
be  viewed  as  a  large  number  of  multiple  reflections  from  the  "granitic- 
basaltic"  layer  interface  in  the  middle  crust;  and  that  the  anelastic  dissipation 
function,  or  Q,  is  frequency  dependent  with  the  anelastic  Q  increasing  with 
increasing  frequency,  and  that  the  low  velocity  zone  absorption  is  srongly 
variable  regionally  (and  also  within  particular  regions)  with  the  low  velocity 
zone  Q  dominating  the  adsorption  along  teleseismic  paths. 
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I.  Introduction 

The  objectiveSof  the  research  being  conducted  are  to:  (1)  Develop 
methods  of  seismogram  synthesis  using  mode  superposition  and  related 
methods,  (2)  Finalize  the  theory  for  source  inversion  by  modal 
decomposition,  (3)  Determine  the  anelastic  characteristics  of  the  medium 
using  known  source  characteristics,  (4)  Interpret  seismic  event 
discrimination  in  terms  of  the  physical  properties  of  the  source, 

(5)  Establish  seismic  event  discrimination  methods  from  formal  inversion 
techniques,  and  (6)  Establish  regional  discrimination  techniques  based 
upon  physical  properties  of  the  source. 

In  this  report  we  describe  research  results  relating  to  seismogram 
synthesis  (item  (1)  above)  and  the  interpretation  of  seismic  event 
discrimination  in  terms  of  the  physical  properties  of  the  source 
(item  (4)  above). 

II.  Seismic  Source  Modeling 

In  order  to  interpret  seismic  event  discrimination  in  terms  of  the 
physical  properties  of  the  source  and  to  be  able  to  establish  new 
discrimination  techniques  we  have  generalized  seismic  source  models 
based  on  relaxation  source  theory  to  include  the  effects  of  non-homogeneous 
initial  prestress  (Stevens,  1980).  In  particular,  we  have  considered 
the  effects  of  strongly  concentrated  prestress  in  the  vicinity  of  the 
shatter  zone  produced  by  an  explosion.  This  work  shows  in  detail  how 
the  effects  of  tectonic  release  can  perturb  the  normal  seismic  radiation 
from  an  explosion.  The  results  of  this  work  are  included  in  the  Appendix  A. 


The  important  results  from  the  work  so  far  completed  are  (1)  The 
spectra  of  P  and  S  waves  radiated  due  to  stress  relaxation  effects  can 
be  strongly  peaked,  with  the  nature  of  the  peaking  being  azimuthally 
dependent  in  general  and  quite  strongly  dependent  on  the  size  and  location 
of  the  initial  stress  concentration,  (2)  The  corners  or  peak  frequencies 
of  the  P  and  S  wave  radiation  are  different  from  one  another  (S  lower) 
and  are  both  shifted  to  higher  frequencies  when  the  stress  concentration 
is  close  to  the  shatter  zone.  In  addition,  the  corner  or  peak  frequency 
value  is  related  to  the  size  of  the  stress  concentration  rather  than  to 
the  size  of  the  shatter  zone,  (3)  the  pattern  of  first  motions  from  the 
tectonic  release,  when  the  prestress  is  inhomogeneous,  is  not  pure 
quadrupole,with  higher  order  multipoles  also  involved.  The  ordinary 
quadrupolar  pattern  predicted  from  a  homogeneous  prestress  can  be  strongly 
distorted  when  the  prestress  is  concentrated  and  can  be  highly  non- 
quadrupole  in  form. 

When  the  prestress  has  an  average  (homogeneous)  component  value  (near 
100  bars  or  so)  and  a  nonhomogeneous  component  corresponding  to  local 
stress  concentrations  in  the  vicinity  of  explosion  produced  shatter  zone 
(with  stress  levels  near  several  hundred  bars)  then,  for  this  "expected 
prestress  environment",  one  can  predict  the  following  consequences  for 
discrimination  and  explosive  yield  estimation: 

1.  Short  period  perturbations  in  the  wave  train  can  be  expected 
to  be  very  complex  due  to  dependence  on  stress  concentration  effects. 
However,  the  perturbations  should  be  small  to  moderate  for  the 
first  cycle  of  the  P  wave  motion,  while  being  significantly  larger 


for  the  later  part  of  the  P  wave  train.  Body  wave  magnitude 
measured  from  the  first  P  wave  cycle  should,  therefore,  be 
minimally  perturbed  by  stress  relaxation  effects.  Further, 
corrections  to  the  first  cycle  of  the  P  wave  for  tectonic  affects 
could  conceivably  be  made  for  purposes  of  yield  estimation. 

2.  Long  period  surface  wave  radiation  is  strongly  perturbed  by 
tectonic  release  effects  within  the  whole  measureable  low  frequency 
band  (i.e.,  from  approximately  5  sec.  to  60  sec.  in  period).  The 
perturbations  in  the  observed  Rayleigh  wave  forms  can  be  such  as 
to  add  to,  or  subtract  from,  the  explosive  generated  Rayleigh 
wave  depending  on  the  orientation  and  magnitude  of  the  prestress 
in  the  vicinity  of  the  explosion.  The  magnitude  of  this  effect 
can  be  very  -(unacceptably)  large.  In  those  cases  where  Love 
waves  are  significant,  so  that  tectonic  release  is  involved,  then 
yield  estimation  using  Mg  can  only  be  made  after  correction  of  the 
Rayleigh  wave  measurement  using  the  observed  Love  wave  to  deduce 
the  size  and  configuration  of  the  tectonic  source.  Such  a 
correction  would  be  much  more  reliable,  when  spherical  shatter  zone 
induced  tectonic  release  is  involved,  than  it  would  be  if  actual 
earthquake  triggering  is  involved. 

III.  Wave  Propagation  Theory:  Synthetic  Seismograms 

In  the  course  of  our  development  of  methods  of  synthesizing  seismo¬ 
grams  in  the  regional  and  teleseismic  distance  ranges  we  have  considered 
two  theoretical  approaches;  in  particular  mode  superposition  using  a 
"locked  mode"  approximation  method  (Harvey,  1980)  and  the  full  wave 


asymptotic  method  (e.g.  Cormier,  1980).  Results  from  the  full  wave 
theory  are  described  in  this  report.  Ti.e  Appendix  B  provides  a  detailed 
discussion. 

The  importance  of  the  full  wave  theory  is  that  it  is  applicable 
at  all  distance  ranges,  incorporates  sphericity,  can  be  applied  to  media 
models  with  velocity  gradients,  is  a  frequency  domain  theory  and  as 
one  consequence  can  incorporate  frequency  dependent  anelastic  absorption 
and  scattering,  and  finally,  it  is  accurate  and  quite  fast  computationally. 

In  our  work  so  far,  the  theory  has  been  fully  developed  and  coded 
for  the  near  and  regional  distance  ranges  from  a  seismic  source.  It 
is  now  possible  to  compute  complete  seismograms  produced  by  rather  complex 
source  models  in  complex  earth  models.  Programs  for  this  purpose  have 
been  used  to  predict  the  radiation  from  explosion  and  earthquake  models 
in  the  near  and  teleseismic  distance  ranges.  Examples  of  these 
computations  in  the  near  regional  distance  range  are  shown  in  the 
Appendix  B,  along  with  the  full  theoretical  development  of  the  method. 

One  application  of  this  capability  has  been  to  study  the  anelastic 
properties  of  the  upper  mantle  with  attention  to  frequency  dependence 
of  the  absorption  (Lundquist,  1980).  In  the  work  to  continue  we  will 
use  this  theoretical  modeling  capability  to  further  study  anelasticity , 
but  also  to  study  both  earthquakes  and  explosions  in  both  the  regional 
and  teleseismic  distance  ranges  to  infer  source  properties  and  to  define 
discrimination  methods. 
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Seismic  Radiation  from  the  Sudden  Creation  of  a  Spherical  Cavity  in  an 
Arbitrarily  Prestressed  Elastic  Medium 

Summary.  We  solve  the  general  problem  of  seismic  radiation  from  the 
sudden  creation  of  a  spherical  cavity  in  an  arbitrarily  prestressed 
elastic  medium.  This  problem  has  direct  application  to  tectonic  release 
due  to  the  creation  of  a  shatter  zone  by  large  underground  explosions. 
In  addition,  however,  the  problem  has  essential  features  in  common  with 
the  more  general  earthquake  source.  Specifically,  it  can  provide  an 
understanding  of  how  an  inhomogeneous  prestress  can  affect  radiated 
seismic  energy  from  a  tectonic  source.  The  problem  is  solved  through 
the  use  of  a  Green's  tensor  integral  equation  in  which  all  quantities 
are  expressed  in  terms  of  vector  spherical  harmonics.  We  obtain  an 
exact  solution  and  an  approximate  solution  to  this  problem.  The 
approximate  solution  may  be  useful  for  the  study  of  other  failure 
geometries.  The  results  agree  with  previous  solutions  for  the  case  of 
pure  shear.  The  case  of  localized  inhomogeneous  prestress  is  examined 
in  detail.  The  primary  result  of  the  inhomogeneous  prestress  is  the 
addition  of  considerable  energy  at  frequencies  above  the  usual  corner 
frequency.  This  causes  an  increase  in  the  corner  frequency,  a  change  in 
the  slope  of  the  spectrum  near  the  corner  frequency  and  in  some  case 
a  strong,  low  frequency,  far  field  spectral  peak.  Peaked  spectra 
always  exist,  however,  near  the  usual  quadrupole  nodes.  Further,  the 
angular  distribution  of  radiation  in  general  will  not  be  pure  quadrupole 
in  nature.  As  an  example  of  a  strongly  inhomogeneous  prestress  case, 
stress  concentrations  near  the  source  due  to  a  point  dislocation  and  a 
point  center  of  compression  are  considered.  The  radiated  waveforms 
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and  spectra  then  vary  greatly  with  angle^  with  the  spectra  peaked 
strongly  at  all  azimuths  in  some  cases.  For  these  cases  the  angular 
distribution  of  first  motions  is  still  relatively  simple,  but  can  be 
deceptive  if  used  to  determine  a  focal  mechanism.  3ecause  of  the 
canonical  relation  of  this  problem  to  the  more  general  earthquake 
failure  problem,  it  is  to  be  expected  that  similar  stress  concentra¬ 
tion  effects  will  occur  for  earthquakes.  Thus,  similar  changes  in  corner 
frequency  due  to  stress  concentration  effects  could  lead  to  errors  in 
earthquake  source  dimension  estimates  and  the  related  existence  of 
spectral  peaks  could  lead  to  errors  in  seismic  moment  measurements. 
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1.  Introduction 

An  earthquake  can  be  considered  to  result  from  the  creation  of  an 
external  boundary  in  a  prestressed  elastic  medium.  Similarly,  at  least 
some  of  the  anomalous  radiation  from  an  explosion  cat  be  considered 
to  result  from  the  creation  of  a  shatter  zone  in  a  prestressed  medium. 

Here  we  consider  the  problem  of  the  sudden  creation  tf  a  stress  free 
surface  in  an  arbitrarily  prestressed  medium  since  this  problem  con¬ 
tains  many  of  the  essential  characteristics  of  the  earthquake  and  explosion 
induced  tectonic  release  problems.  In  particular  it  shows  clearly 
and  simply  what  effects  inhomogeneous  prestress  will  have  on  the 
observed  radiation  field  We  will  show  that^problea  nay  be  treated  as 
an  initial  value  problem  or,  equivalently,  as  a  stress  pulse  problem. 

We  solve  this  problem  exactly  for  the  case  of  a  spherical  cavity 
using  vector  spherical  harmonics  in  the  elastic  Green’s  tensor  integral 
equation  in  the  frequency  domain.  The  result  is  a  relatively  simple 
solution  for  the  radiation  field  involving  nothing  mere  than  the  inverse 
of  a  number  of  two  by  two  matrices  and  an  inverse  Fourier  transform. 

The  solution  is  valid  for  any  arbitrary  homogeneous  cr  inhomogeneous 
prestress  and  is  valid  in  the  near  field  or  far  field.  Previous  solutions 
to  spherical  source  problems  have  all  been  obtained  f:r  uniform  stress 
fields.  Hirasawa  and  Sato  (1963)  found  an  exact  solution  for  the  sudden 
creation  of  a  spherical  cavity  in  a  pure  shear  stress  field.  Randall 
(1964a)  used  a  spherical  inclusion  in  which  the  material  inside  suddenly 
underwent  a  phase  change  to  become  the  same  as  the  external  material  to 
model  deep  earthquakes.  He  applied  this  to  a  uniform  compressive 
(1964b)  and  a  uniform  shear  (1966)  prestress.  Archamteau  (1968)  used 


4. 


a  growing,  propagating,  transparent  spherical  cavity  in  a  uniform  shear 
field  as  a  model  for  earthquakes.  Minster (1973)  refined  this  model. 
Burridge  and  Alterman  (1972)  found  an  exact  solution  fcr  a  uniformly 
growing  spherical  cavity  in  a  pure  shear  field.  Burridge  (1975)  used 
this  solution  as  a  test  for  Archambeau's  (1965,  1972;  transparent  source 
approximation.  Koyama  et  al.  (1973)  found  an  exact  solution  for  the  sudden 
creation  of  a  fluid  filled  cavity  in  a  pure  shear  field.  Minster  and 
Suteau  (1977)  examined  the  relationship  between  a  grcving  spherical 
source  and  a  growing  circular  dislocation  with  the  seme  growth  history. 

The  case  of  inhomogeneous  prestress  is  particularly  interesting 

since  almost  all  previous  source  models  have  used  the  boundary  condition 

of  homogeneous  shear  or  its  equivalent.  Actual  stress  fields  in  the 

earth  are  likely  to  be  highly  inhomogeneous.  Arckambeau  (1968,  1972) 

attempted  to  allow  for  a  localized  prestress  by  solving  the  static 

initial  value  problem  for  a  homogeneous  shear  field,  but  then  limiting 

the  effective  source  region  to  be  within  a  radius  Rc  for  his  dynamic 

calculations.  Snoke  (1976)  showed  that  there  were  problems  with  this 

approximation.  Bache  and  Barker  (1978)  attempted  to  determine  stress 

variations  in  the  earth  during  the  1971  San  Fernando  earthquake  by 

■?  v  '  e 

using  a  solution  for  a  transparent  growing  sphere  in  a  -ghoar  shear 
field,  but  allowing  the  magnitude  of  the  stress  field  to  vary  during 
growth.  We  find,  in  an  exact  calculation,  that  the  maim  effect  of  the 
inhomogeneous  prestress  is  appearance  of  considerable  seismic  energy  at 
frequencies  greater ‘than  the  usual  quadrupole  corner  frequency.  A 
slightly  inhomogeneous  prestress  results  in  a  nearly  quadrupole  dis¬ 
tribution  of  radiation,  but  with  anomalous  radiation  in  and  near  the  nodes. 
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A  stress  concentration  close  to  the  source  results  in  a  dramatic  increase 
in  corner  frequency  and  the  appearance  of  pronounced  spectral  peaks  at 
all  points  of  observations.  The  waveforms  produced  by  a  concentrated 
prestress  are  substantially  different  from  those  produced  by  a  uniform 
prestress.  The  angular  distribution  of  first  notions,  however,  nay  or 
may  not  appear  to  have  a  quadrupole  distribution. 

This  paper  is  the  first  attempt  to  explicitly  calculate  the  effect 
of  an  inhomogeneous  prestress  on  the  spectra  and  waveforms  of  a  dynamic 
seismic  source.  The  spherical  cavity  problem  can  be  solved  exactly  and 
is  the  canonical  problem  which  can  be  used  to  determine  some  of  the 
features  of  failure  in  other  geometries.  The  results  should  be  important 
to  observational  seismologists  since  the  corner  frequency  is  often  used 
‘to  estimate  the  size  of  a  seismic  source  and  the  results  of  this  study 
indicate  large  variations  in  the  corner  frequency  as  a  function  of  local 
prestress  inhomogeneity,  location,  and  magnitude.  In  some  cases,  for 
example,  increase  in  the  apparent  corner  frequency  occurs  and  this  can 
cause  a  gross  underestimation  of  source  dimensions.  Further,  the 
complexity  of  the  spectrum  that  can  be  produced  can  make  the  corner 
difficult  to  define  observationally  and  peaking  of  the  spectrum  can 
lead  to  uncertain  or  erroneous  measurements  of  moment.  On  the  other 
hand,  the  theory  given  here  should  prove  useful  for  the  study  of  spatial 
variations  of  stress  in  the  earth. 

We  are  simplifying  the  seismic  source  problem  in  this  paper 
through  the  use  of  instantaneous  creation  of  the  failure  surface  and 
through  the  neglect  of  any  consideration  of  the  material  inside  the 
cavity.  In  this  way  we  can  examine  the  effect  of  an  inhomogeneous 


stress  field  on  seismic  radiation  independently  of  the  complexities 
of  finite  rupture  velocity  and  certain  boundary  effects.  The  formidable 
general  problem  of  failure  in  a  prestressed  medium  in  which  the  growth 
rate  is  determined  by  the  prestress  and  the  problem  is  coupled  at 
the  boundary  to  the  material  inside  the  cavity  has  been  discussed  in 
detail  by  Archambeau  and  Minster  (1978). 
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2.  Basic  Relations 

The  problem  of  the  radiation  released  when  a  cavity  forms  in  a  prestressed 
medium  is  conveniently  expressed  in  terms  of  a  Green’s  tensor  integral 
equation.  Consider  first  the  cavity  of  arbitrary  shape  shown  in  figure  1. 

We  assume  that  the  medium  is  prestressed  and  that  the  cavity  forms  at  time 
t  “  0.  The  medium  is  then  displaced  from  equilibrium  initially  and  relaxes 
to  a  new  equilibrium  state.  The  displacement  field  satisfies  the  elastic 
equations: 


^•u(x,t)  -p—  u(x,t)  -  0 
3t 


C2.1) 


3  3 

where  L- u(x,t)  =  -r —  —  u. 

- -  3x^  ijk£,  3x^  2, 


(2.2) 


is  the  elastic  operator,  is  the  constitutive  tensor  for  the  medium 

external  to  the  cavity. 

The  Green’s  tensor  satisfies  the  equation: 


a 

L*G(x,  x  ,t,t  )  -  p  — 7r  G(x,x  ,t,t  )  ■  -  I  S(x-x  )  SCt-t  ) 
~o’  ’  o'  gt2  ~  ~  ~o’  *  o  ss  ~~o  o 


(2.3) 


We  will  solve  this  problem  in  the  frequency  domain  with  the  aid  of 
vector  harmonics.  We  define  the  Fourier  transform  and  its  inverse: 


u(w) 


/•: 


(t)  e'1”'  dt 


u(t) 


27T  J  “ 


e  \  i“t  , 
(u>)  e  dw 


(2.4) 


The  transformed  displacement  field  satisfies 


L*u(x,uj)  +  pw  u(x,oj)  -  0 


The  transformed  Green's  tensor  satisfies: 


(2.5) 


i.*£(x,x  ,ai)  +  p<i>  G(x,  x  ,w) 


-  I  6(x-x  ) 
z  ~  ~o 


(2.6) 


It  will  prove  simpler  to  use  vector  notation  rather  than  indicial 
notation  throughout  this  paper,  and  we  define  the  stress  operators: 


Cijk£  3xfc  U£  ~  Tij 


3x^  ^Jtm 


(2.7) 


(2.8) 


The  Green's  tensor  integral  equation  for  the  general  case  of  growing 
phase  boundaries  has  been  derived  by  Archambeau  and  Minster  (1978) ,  For 
the  case  of  the  instantaneous  creation  of  a  boundary  it  reduces  to: 


u(x,t)  =  fdto  J 


[G-T(u) 


-  u*  T(G)  ]  *ndA 

V  $5  »  C 


■m 


(2.9) 
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The  last  term  contains  the  effect  of  the  initial  conditions.  The  surface 
integral  contains  the  response  of  the  surface.  The  first  term  in  this  integral 
will  vanish  for  a  stress  free  surface. 

The  frequency  domain  integral  equation  is: 


u(x,a>)  -  -J  [u-T(G)  -  G-T(u)]* 

+  iaj  /  pu  'G  dV 
I  ~  »  o 
JV 


(2.10) 


where  u*  =*  u(t**0),  the  initial  displacement  from  equilibrium.  Again  the  term 

containing  T(u)  vanishes  for  a  stress  free  surface. 

»  * 

The  initial  value  term  in  equations  (2.9)  and  (2.10)  can  be  quite  difficult 

to  evaluate.  The  problem  can  be  made  much  simpler  by  defining  the  "stress 

pulse  equivalent"  to  the  initial  value  problem.  We  use  the  fact  that  since 

u^*  is  the  difference  between  two  static  equilibrium  fields  (the  prestressed 

state  without  the  cavity  and  the  final  equilibrium  state) 

L- u*  =  0  (2.11) 

&  ~ 

Multiplying  this  equation  by  G  and  multiplying  equation  (2.6)  by  u*, 

fz/ 

subtracting  and  integrating  overall  space  external  to  the  cavity  and 
substituting  in  equation  (2'.‘10)  we  find  that  for  a  stress  free  cavity  of 
arbitrary  shape: 


u'«T(G)*n 

/V  21  S 


dA 

o 


(2.12) 


where  u’  is  the  relative  displacement  field 


u’(x.iu)  »  u(x,u)  -  U*(x)/i(D 

T(u*)  .  n  is  the  traction  drop  on  the  cavity  surface  between 
the  initial  and  final  state.  Since  the  final  tractions  are  zero  for 
the  stress  free  cavity,  T(u*)  .  n  is  just  the  initial  tractions  on 
the  surface  due  to  the  existing  prestress.  Comparing  equation(2.12) 
with  equation  (2.10),  we  see  that  the  initial  value  problem  is  equivalent 
to  a  stress  pulse  equal  in  magnitude  to  minus  the  initial  tractions 
applied  to  the  cavity  surface  at  time  t  =  0  .  While  the  relaxation 
problem  is  perhaps  more  naturally  expressed  as  an  initial  value  problem, 
the  solution  is  much  easier  when  expressed  as  a  stress  pulse  equivalent. 
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The  infinite  space  Green’s  tensor  in  the  frequency  domain  is  given 
by  Ben-Menahem  &  Singh  (1968)  in  terms  of  vector  spherica.1  harmonics: 


G(x,  x  ).*  -  -o— ^ ^ 
&  ~  ~o  4iru  .  -- 


*-1.1.0 


22+1  ^  (2-r) 

*(*+!>  (£+n) 


Jn.-A  ?  (r  )  M~C  (r) 
(2+n) !  »v/.3  o  /vs,mv 


+lk(ro>  ^(r>  +  i(t+i)  (!)3-' 


(2.13) 


where  e  =  sgn  (r-rQ)  and  r  =|x|  . 


V  x  (e  *-)  ;  N"  = 


*r*p,L- 


1  ,± 

k  Va 


k  -  -  •  k  *  ^ 

a  a  •  kB  B 


Here  a«B  are  the  P,  S  velocities  of  the  medium 


£  =  Y2m(6»^ 


(  ^(kctr)  1 
)  hn(2)(kr)) 


While  is  similarly  expressed  using  kg  . 

The  operation  *  refers  to  complex  conjugation  of  the  angular  part 
of  the  function,  while  the  Hankel  function  remains  unchanged.  Explicit 
forms  for  M  ,  N  and  L  are  given  later  (see  equation  .3.9)  -  (3.13)). 
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3.  Solution  for  Arbitrary  Prestress 


By  expanding  all  relevant  quantities  in  terms  of  vector  spherical 
harmonics  and  using  the  convenient  orthogonality  properties  of  these 
functions,  we  can  find  an  exact  solution  for  the  radiation  field  for  the 
sudden  creation  of  a  spherical  cavity  in  an  arbitrarily  prestressed  medium. 
We  need  to  solve  the  equation: 


u  *  /  u-T(G; 

Arf  /  -V  A/  V 


T(G)  *n  dA  + 


(3-D 


Here  n  is  the  outward  normal  e^  at  the  spherical  surface  . 
u1  may  be  either  the  initial  value  term: 


I 
u  - 


+  iup 


f 


G  dV 


(3.2) 


where  V  is  the  volume  external  to  the  sphere,  or  the  stress  pulse  term: 


i1  =  ^-  /  G-T(u*)’S  dA 
iwj  &  zs  ~ 


(3.3) 


The  stress  pulse  term  is  easier  to  evaluate.  We  will  now  solve  equations 
(3.1)  and  (3.3)  for  a  general  stress  field. 

Any  vector  function  can  be  expanded  in  the  following  way: 


00  i 


+  a.^W.C. 

x.n  i 


(3- A) 
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P.  ,  B.  ,  C.  are  vectors  defined  by  (see  Morse  and  Feshbach  (1953)  or 
-s.£m  ~£m  'v£m 

Ben-Menahem  and  Singh  (1968)): 


>  £m 


^r  Y£m<9’*> 


(3.5) 


/{,  ( £/fl7  m  -  £eQ  3Q  +  e^  gin  g  3(f)  J 

/£(£+!)  C£m  =  £  eQ  sln  0  "  e(j)  Jq  j 


Y£n<0 »♦>  =  PJLmCcos  6)  eimt 


where  P^m(x)  *s  the  associated  Legendre  function.  These  vectors  have  the 
following  orthogonality  relations: 


/  P.  -B.  ,  ,  dft  =  f  P„  -C.  ,  tdQ  =  f  B.  -C0  ,  ,dft  =  0 
J  ~  £m  ~£  m'  J  ~ £m  -Jl'in  J  -*£m  ~£  m 


and 


(3.6) 


/P„  "P-,  ,dn  =  f  B„  •  B  ,  ,dfi  -  /*  C„  Cn>  ,dft  =  6„„,<$  ,«r 

^  ,vi.m 'vi '  m'  J  *»£,m--£'ni'  ££'  mm'  £m 


n 


(£+m) ! 
£m  ~  2£+l  (£-m) ! 


Either  by  inspection  or  by  using  the  orthogonality  relations  we  can 
expand  u*  in  the  form: 
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u>,r,e,*)  =  T  E  Wn(r,«)JPlm(e,«) 

£=0  np-r 

+  d?  (r,u)  B.  (6,4.)  +  d?  (r,w)  C.  (6,6)1 
in  £  m  xin  ~  ju  z.  j 

We  can  also  expand  the  unknown  displacement  field  u  in  the  form 

u(«,r,e,«  -  E 

£=0  m=-£ 


(3.7) 


(3.8) 


+  bL(r*“)  Jan<0'*>  +  bL(r»“}  -W9’*5} 


The  Green's  tensor  for  this  problem  is  given  by  equation (2 .13) .  In 
terms  of  the  vectors  P,  B,  C  we  can  write  (Ben-Menahec  and  Singh  (1968)) 


Jim  81^  ~£m 

iSm=  m+  *3(y>  ~£m 

~£m  g4(x)^£m+  4(x)  Jim 


(3.  S) 


where  y  =  k^r,  x  =  k^r.  The  vector  M  represents  toroidal  waves.  The 
vectors  L  and  represent  P  and  S  spheroidal  waves  respectively. 

We  also  need  the  associated  tractions 

iq£m).S  -  „!>*<*>  clm 

T<!>«  -  »t J 

l(!+1>  («  f  £< ‘  »(h4<’0  *1*  +  ^5 (x>  i  J 


(3.10) 
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where 

g^(y)  =  /Ta+i)  f*(y) 

87*  (y)  =  £(£+1)  -^1 

i.  y 

g£±(y)  "  /£<£+D  ^fjf(y)  + 

4+  —  * 

g£  (x)  =  f'  (x) 

5+  _ _ ff  (*) 

g.  (x)  =  vCu+1)  — - 

X 

hf  (y)  =  k3  /mm  y  F^(y) 
hf  (y)  -  2kg  £(£+1)  F^(y) 

h£  ^  =  k6  ^  F£3^  (3.11) 

h£  (x)  =  2ka  F£A(x)  £(£+1)^f  j 
h^(x)  =  2ka  /£  (£+1)  F^(x)  £(£+l)(f  y 

F£1(Z)  "  (2£~i)  (2£+l)  (2£+3)  [  <£-D  <2£+3)  f*_2  (z)  -  C21+1)  fj(z)  -  (£+2)  (2£-l)  f*+2U) 

F£3^z)  "  (2£-l)  (2£+l)  (2£+3)  [  2^£  _1)  (2£+3)  f^_2(z)  *  3(2£+l)f£(z) 

+  2£(£+2)(2£-l)  f£+2(z)  j 
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Fu(*>  -  W-tttt+inaa)  [»<*-»  (21+3)  4_2M  *  <£+l)tt+2)(2£-l)  42(z) 


-  (2£+l)  (  2£2  +  2£-  1  +  — ~f-8 

{  26^ 


(2£-l)(2£+3) 


where  f~  are  spherical  Bessel  (+)  and  outgoing  Hankel  (-)  functions. 

(There  is  a  misprint  in  Ben-Menahem  and  Singh's  definition  of  which 
is  corrected  here.) 

The  surface  integral  in  equation  (3.1)  can  now  be  evaluated.  Because 

of  the  orthogonality  relations  (3.6)  the  surface  integrals  give  either 
2 

R  $J  or  zero.  The  integral  equation  has  now  been  reduced  to  a  set  of 
algebraic  equations.  In  fact,  the  individual  terms  in  the  eigenfunction 
expansion  are  almost  completely  uncoupled.  The  remaining  equations  for  each 
value  of  £  are: 


(r) 


-ik/ 

£(£+1) 


(s£~(r)  b£m(R)  h>  +  8*4’(r)  b^(R>  ^+(R) 
+  4’(r)  bL(R)  h>  +  8i"(r)  btm(R)  h*+ 


(R) 


(r) 


bL<r)  -  bL(R)  h£+(R)  +  b£m(R)  h£4+(R) 


(3.12) 


+  b£m(R)  hf(R)  +  8^'(r)  b^(R)  ^+(R)J 


+dL(r) 


b£m(r) 


imr  lhM  hL™  h£1+(R) 


+  d£m(r) 


(The  arguments  of  the  g  and  h  functions  are  understood  to  be 


multiplied  by  k  or  kQ  as  defined  in  equation  (3.11)). 
a  p 

We  can  now  evaluate  these  equations  at  r  =  R  (the  cavity  radius) 


and  have  the  solution  to  the  problem. 


The  coefficients  of  are  uncoupled  from  the  others,  so  we  have 
-ik^R 


immediately  letting  — ■  Q$ 


bL(R> 


1  -  Q£  gJ"(R)  hJ+(R) 


(3. 


and  at  any  point  r 


23 


X  2 

The  coefficients  b*m  and  b*  are  coupled  for  £  f  0.  The  solution  for  these  is 


'bLm 


»Lw 


(3.15) 


where  Au  -  g“'(R)  h*+(R)  +  g*~(R)  h*+(R) 


A12  ■  g£"00  hf'OO  +  g*~(R)  h^OO 


(3.16) 


and 


A2i  -  gf“(R)  hf'(R)  +  g£~(R)  h*+(R) 


A22  =  4  (R)  hf<R)  +  %~(R)  hf(R) 


bJ-(RA 

£m 


>bLR)y 


-l 


hr 


(3.17) 


This  requires  the  inversion  of  only  a  2  x  2  matrix.  Substitution  into 
equation  (3.12)  gives  the  exact  solution  to  this  problem. 

The  case  £= 0  must  be  treated  separately  since  3  ■  C  =0.  Setting 

r  —00  -'OO 

2  1 

b  ■  0  we  can  immediately  solve  for  b 
oc  oo 
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bL(R) 

oo 


dL(R) 

oo 


l-qo  840'<H  h4w 


(3.18) 


bnn(r) 

OO 


Qog*"(r)  h*+(R)  d*  (R) 

d1  (r)  +  _°_2 - 2 - 22 - 

oo  ,  , . 

1  -  Q  g*(R)  h^(R) 
o  o  o 


(3.19) 


I 


/ 

The  factors  of  £,(£,+1)  in  Q  and  h  cancel. 

O  O 

Inspection  of  the  solution  (equation  3.12)  shews  that  the  displacement 
field  consists  only  of  outgoing  spherical  waves.  If  desired,  the  solution 
can  be  written: 


u(x,u) 


£ 

£.=0 


m=-£. 


a.  (u>)  M  + 
£.m  'vim 


L. 


*s.m 


where  the  coefficients  a  ,  B.  ,  y  are  determined  by  comparison  with 

Jem  Jem  Jem 

equation  (3.12).  The  first  term  is  a  toroidal  wave,  the  second  a  spheroidal 
shear  wave  and  the  third  a  spheroidal  compressionel  ware. 
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4.  Pure  Shear 

The  case  of  the  creation  of  a  spherical  failure  surface  in  a  pure 
shear  field  has  been  used  as  an  element  in  a  model  fcr  earthquakes  by 
Archambeau  (1968)  and  Randall  (1966).  Archambeau  (IS 72)  used  this  as  a  model 
to  predict  anomalous  stress  wave  radiation  from  explosions.  Randall's  solution 
is  equivalent  to  the  stress  pulse  solution  which  is  used  here  as  the  Initial 
value  term.  This  problem  was  solved  exactly  by  Eirasawa  and  Sato  (1963). 

Their  results  are  reproduced  in  Koyama  et  al.  (1973).  K^rama  et  al.  solved 
the  problem  of  the  creation  of  a  fluid  filled  cavity  in  a  pure  shear  field. 

They  report  little  difference  in  the  waveforms  for  the  tvo  cases. 


It  is  possible  to  do  this  problem  in  two  ways.  An  expression  for  u*  is 

given  by  Landau  and  Lifschitz  (1970.  p.  24).  After  some  effort  it  can  be  shown 

that  this  can  be  written  in  terms  of  vector  spherical  harmonics  with  only  the 

coefficients  of  B„  and  P.  non-zero.  The  initial  value  integral  (equation 
'vam  am 

3.2)  can  then  be  performed.  This  leads  to  a  very  messy  expression.  It  is  much 
easier  to  use  the  stress  pulse  solution  and  evaluate  equation  (3.3). 

T(u*)*h  is  easily  found.  The  initial  normal  tractions  are  simply  those  of  a 

A/  ^ 

pure  shear  stress  field  resolved  onto  a  spherical  surface.  The  final  normal 

tractions  are  zero.  So  T(u*)*n  =»  a .n.  with  c,,  +  c„_  +  r_,  =  0  for  the  most 

m  ~  ~  ij  j  II  j3 

general  pure  shear  field.  We  need  to  write  these  it  tarns  of  vector  spherical 
harmonics.  After  some  algebra  it  can  be  shown  that  the  m:st  general  pure 


shear  field  can  be  written: 
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"  Jt2  A2*  fcm  +  ’I i2m) 


(4.1) 


where  are  coefficients  determined  by  the  choice  of  cartesian  stress  tensor. 


In  terms  of  the  cartesian  components  of  the  shear  stress: 


A22  "  12  (all  "  a22)  +  6i  a12 

A21  “  I  c13  +  3i  °23  . 

A20  “  a33 

A2-l  "  “  2ct13  +  i  °23 
A2-2  ”  2(°ll  "  a22}  '  I  °12 


(4.2) 


Now  we  can  evaluate  equation  (3.3) 


u1  =*  ~  /  G-T(u*).n  dA 
iio  /  « &  - 

r* 


The  non-vanishing  part  of  G  is 


~  ,  ~r — —  y;  j2~-mj-j-  (Ntcx)N"6c)  +  6(-)  l^6c)l'(^; 

«’/«o  4ttu  6  u~,2  (2+m) !  y~«2m  — *o  —2m  ~  \  a  /  '2m  — o  —2m  — ' 


G(x  ,x  ) 


(4.3) 
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Evaluating  the  surface  integral  we  find 

2 

-1-  t£^[4i2.+4j2.] 


where,  using  the  notation  of  equation  (3.11) 


■4.  ’  '  fie  a2«  [^2+(yo> +  4  8f  (r Sj  h  (y) 

+  6(f  )  +  4  82+(v)  824"<x)] 

(4.5) 

d2m  ‘  -  STB  A2m  [^2+(y<,)  +  4  H  (yo>)  4' <y) 

+  6  (f  )  («2+<V  +  4  4X>)  «T« 

x  =  k  R  x  =  k  r  y  =  k„R  y  =  k„r 

o  a  a  Jo  &  s 


These  coefficients  can  be  substituted  into  ecuaticn  (3.17)  to  obtain  the 
exact  solution  to  this  problem.  In  the  far -field  approximation  g2  and  g,. 
vanish.  The  solution  then  corresponds  to  a  radial  P  wave  and  a  transverse 
S  wave  as  it  should.  The  far  field  P  and  S  spectra  for  the  initial  value/ 
stress  pulse  (eq.  4. A)  solution  and  the  exact  solution  are  shown  in 
figure  2.  The  initial  value  spectrum  contains  a  large  cumber  of  dips  not 
present  in  the  exact  solution.  Burridge  (1975)  pointed  out  that  these  dips 
result  from  the  stress  discontinuity  at  the  enc  of  the  pulse.  This  is  a 


consequence  of  the  fact  that  without  the  inclusion  of  the  surface  term 
(eq.  3.1)  there  are  in  effect  no  boundaries  in  the  medium.  All  disturbances 
are  then  propagated  at  the  corresponding  velocity  of  the  medium.  The  resulting 
pulses  will  have  a  duration  equal  to  the  travel  time  across  the  cavity.  The 
farfield  P  and  S  waveforms  are  shown  in  figure  3  for  both  the  initial  value 
and  exact  solutions.  The  solution  agrees  exactly  with  the  solution  of 
Hlrasawa  and  Sato  (1963).  Their  spectra  and  waveforms  are  reproduced  in 
Koyama  et  al.  (1973)  and  match  figures  2  and  3. 
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5.  Radiation  from  Non-Uniform  Stress  Fields 

We  want  to  consider  the  radiation  from  a  stress  free  sphere 
spontaneously  created  in  an  arbitrarily  prestressed  elastic  medium. 

This  should  be  a  good  model  for  predicting  the  tectcnic  release  from 
explosions.  It  should  also  be  a  good  model  for  the  study  of  spontaneous 
failure  in  an  inhomogeneous  stress  field.  Minster  and  Suteau  (1977) 
have  examined  the  question  of  failure  geometry  in  sene  detail.  The 
main  differences  between  the  planar  source  and  the  spherical  source  is 
that  for  the  planar  source  the  amplitude  is  reduced  by  a  factor  of  3 
and  the  waveform  is  smeared  out  over  a  longer  period  of  time  due  to 
interference  of  spherical  waves  leaving  different  parts  of  the  surface 
at  different  times.  This  smearing  effect  can  also  be  seen  in  the  numerical 
results  of  Madariaga  (1976).  While  we  do  not  expect  exact  agreement 
between  the  radiation  produced  by  an  earthquake  and  radiation  produced 
by  a  spherical  cavity,  the  cavity  can  provide  information  on  general 
features  such  as  the  angular  distribution  and  approximate  waveforms. 

Since  the  spherical  cavity  problem  is  exactly  solvable,  it  provides  a 
useful  tool  for  the  investigation  of  the  general  problem  of  failure 
in  an  inhomogeneous  stress  field. 

The  first  question  to  be  answered  is  the  nature  of  different  kinds 

of  prestress  which  are  physically  acceptable.  Any  region  of  the  earth, 

especially  seismically  active  regions  will  contain  many  dislocations. 

These  dislocations  will  act  like  some  complicated  distribution  of 

(equivalent)  static. body  forces  in  the  medium.  The  initial  field 

satisfies  L  .  u*  »  f  and  the  final  field  satisfies  L  .  u^  =  f  ,  so 
+**  '■u  *** 

the  difference  field  satisfies  L,  .  u*  =  0  ,  the  static  elastic  equation 


(  — 
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with  no  body  forces  acting.  This  difference  field  must  be  finite  at 
infinity.  Sicne  we  are  interested  in  the  region  outside  r  =  R  ,  the 
solution  representation  need  not  be  finite  at  r  =  0  .  The  cost  general 
static  displacement  field  satisfying  these  conditions  can  be  expressed 
in  terms  of  the  following  eigenvectors  (Ben-Menehem  and  Singh,  1968): 

A.  ■  r~*"2  I'*  (»+!>,»,„  -  CM-U.P.J 

Jim  "  r"‘  t  1  <V_t)  Jim  +  +i  - 

Jim  ■  r'*'1  Ji„ 


where 

,2 

Y  -  4(l-a)  -  -/-g-  ■ 
a  -  e 

The  displacement  field  is  then 

“‘<r>  ‘  S  +  ‘u  •Eim  +  £i»  Jim  «•« 

are  now  in  a  position  to  compute  the  radiation  field,  me  easiest 
way  to  do  this  is,  again,  to  use  the  stress  pulse  solution 

=  J~  •  _T(G)  •  ndA  +  _uj 

l 


u 


(5.3) 
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where 


T  1  r  *  r. 

u  =  -r-  /  G  •  T(»  )  •  adA 


(5.4) 


*  /\  * 

We  need  to  comoute  T(u  )  *  e  for  u  give 

*  -*w*  J-  /V 

stress  field  is  given  in  spherical  coordinates  by: 


n  in  equation  5.2.  The 


.  3ur 

Trr  =  Xe  +  2lJ  IT 


3u6  1  3ur 

T9r  =  U[3r~  +  r(39  u8)] 


3ui>  1  3u 

hr  ‘  +  73ss  - sin-  V: 


(5.5) 


where 


e  *  sT1  +  7  (3^-  +  2  ur)  +  r"H7e  is*  +  "7^-  ue 


We  consider  any  function  of  the  form: 


u(r,6,$)  =  E[£l«n(r>^A,n(9**>  +  f2£m(r)  im(?  >f)  +  f31m(r)  £ln(-  ’*}] 

JL,m 


(5.6) 


Using  the  expressions  for  _P  ,  _B_  and  jC,  (equations  3.5)  and  the  equation 
for  Y£jr(6,^)  : 


r  -\2 


1 36' 


+  cot?  -Ka  +  - T~  +  ^(^+D 

^  sin  9  3<? 


T£_(f,c)  *  C 


(5.7) 
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we  have  after  some  algebra: 


•  w3fUm  2fUm 


f2£ms  ,  3fUm, 


t(u)  •  er  =  £  W-ar2  +  — *  - 


+  y(^|^  +  /&U+TT  Bia(e,«) 


.  ,9f 3£m  .  f3£m,  „  /A 

+  p(-37 - -) 


(5.8) 


comparing  with  expression  (5.2)  we  find 


*<*.*>  '  *t  ‘  E-ESm(e-«  K,  r"*"3  '2''<,1+1>  <i+2>l 

n  _  ' 


,+  ^!Sr-X'-1  [(y+Ji-i)  (\ Z-212  (A+m))  +  AA(A+1)(2A-1)]} 


+  /A(A+1)  8^(6, <|.)  |aim  r  £_3  2u(-(£+2)) 
+  -iS  r"^1  y(2i2-y)| 

+  A  (1+1)  0^(0, 4>)  |fin  r-£-2  u(-(A+2)j. 


(5.9) 


Comparing  this  with  results  from  the  last  section,  we  find  that  pure  shear 


corresponds  to 


a_  =  4 
2m  5  2m 


c2m  ‘  *  k2» 


(5.10) 
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where  are  constants  depending  on  the  stress  iiagultuie  and  orientation. 

Using  the  Green's  tensor  (eq.  2.13)  for  r  >  with  tq  *  R,  and  using 
the  notation  of  equations  {"3.13),  we  can  now  evaluate  equation  (5.4)  for  the 
most  general  initial  stress  field,  we  have: 


T  _2  “  ,  £ 

I  -R  V"  1  v" 

-  “  £=o  ^£+1)  ±Li 


(5.11) 


f£m  8 


£~(k8r)[S£l 

r)  [S£2(R)  8 


(R)  si  (k„R)  C4m(0,i) 


+  a£m  h  (k8r)  S£2 (R)  8i  (kSR)  +  Sp^R>  8c  <k=R> 


£4  s£ 


+  C£m  8£_(k2r)  S£3(R)  8£  (k6R)  +  S£5(R) 


+  a£m  8£  ^k8r^  S£2(R)  g£  (kS 


+  C£m  8£  ^k8r^  S£3^  8£  ^kgR^  +  SJ5^R^  8 


R)  +  s15(R)  gf+(k£R)  ?la'e,6) 

R)  +  s£4(R)  g?* (k.R)  3£a;T,o) 

— 

R)  +  s£5(R)  g*  (Jc,R;  3ic  :rO 


+  a£m  8£  (V}  1 (*+1) 


+  c£m  g£  (kar)  *(£+1) 


+  a£m  g£  (kar)  £ (£+1) 


+  c£m  8£  (kar)  £<*+1) 


(- )’  [ 
(*>■[■ 


s£2(R)  8£  (koR)  +  5 


.,(7.)  r'  Ck3R)Jp4B(ei< 


S£3(r)  g£  (kaR)  +  s. . (2.)  rr  (M)  jtB(e.« 


s£2(R)  g£+(kQR)  +  s£4(S)  rT  (k^R) 


B.  (0,< 

i  -vj.nl 


s£3(R)  g£  (kaR)  +  s£5(?0  5I  (k 


(6,*) 
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where  the  static  s  functions  are  given  by: 


sn  ' 


2y(£+l)(£+2) 

S£2  =  r£+1 

_  (y+£-1)  (A£-2£2(l+y))  +  A£(£+l)  (2£-l) 


_  -  2y( £+2 )/£(£+!) 
S£4  £+1 

K 

y (2£2-y) /£(£+l) 

s£5  £-1 

Y  R 


This  is  the  initial  value  field  corresponding  to  a  stress  pulse  applied 

at  r  =  R  at  time  t=0  in  an  infinite  medium.  The  first  term  in  this  expression 

corresponds  to  a  toroidal  shear  wave.  The  next  four  terms  corresponds  to 

a  spheroidal  shear  wave  and  the  last  four  terms  correspond  to  a  congressional 

2-  5- 

wave.  In  the  farfield  =  g^  =0.  Then  the  remaining  terms  correspond 

to  a  pure  radial  P-wave  and  a  pure  transverse  S  wave.  The  coefficients  of 

12  3 

P,  B  and  C  correspond  to  d.  ,  d  „  and  d.  in  equations  (3.14)  and  (3.17). 

^  ^  'V'  Jem  Jem  Jcia 

Thus,  (5.11)  determines  the  complete  solution  for  the  radiation  field  due  to 
creation  of  a  spherical  cavity  with,  the  most  general  prestress. 

The  coefficients  a.  ,  c.  ,  and  f„  are  independent  until  a  specific 

Jem  Jem  Jem 

presf-ess  is  defined.  These  coefficients  depend  on  the  static  field  only 

and  are  independent  of  frequency.  There  will  therefore  be  three  independent 

waves  -  a  toroidal  wave  and  two  spheroidal  waves  for  each  value  of  £.  The 

spectra  for  £  =  2  are  shown  in  figure  4.  The  two  spheroidal  waves  corresponding 

to  coefficients  a„  ,  c.  are  labeled  PI,  SV1,  and  P2,  SV2,  respectively.  The 
£m  £m 
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toroidal  mode  is  labeled  SH.  It  is  interesting  to  note  that  the  £  -  2 
modes  consist  partially  of  peaked  spectra  and  partly  of  spectra  which  are 
flat  in  the  low  frequency  limit.  The  linear  combination  of  these  (eq.  5.10) 
that  corresponds  to  pure  shear  has  a  flat  spectrum. 

The  £  =  2  waveforms  are  shown  in  Figure  5.  These  look  much  different 
than  the  waveforms  for  pure  shear,  but  in  fact  the  linear  combination  of  these 
given  in  equation  (5.10)  is  the  pure  shear  pulse. 

The  waveforms  for  higher  £  values  have  more  and  more  oscillations. 

An  example  is  given  in  Figure  6  for  the  £  =>  4  waveforms.  The  result  of  this 
is  that  the  spectra  become  more  and  more  peaked  for  higher  order  £  values 
and  the  maximum  frequency  increases  slightly  for  each  increment  in  £.  The 
slope  of  each  peak  as  a  +  0  increases  with  £.  The  slope  of  each  spectrum  as 
ui  ■+  ®  is  -2  for  all  £.  Figure  7  shows  the  individual  PI  spectra  for  £  =  3, 

4,  5  and  6.  The  magnitude  of  these  peaks  will  vary  with  angle  as  the 
corresponding  vector  spherical  harmonic.  It  will  therefore  vary  more  rapidly 
with  angle  than  pure  shear. 

These  results  shed  some  light  on  the  old  argument  about  low  frequency 
spectral  peaks  in  relaxation  problems.  Archambeau  (1968)  solved  the  problem 
for  the  radiation  from  a  spherical  cavity  in  pure  shear,  but  he  made  an 
approximation  for  a  more  localized  stress  field  by  truncating  the  volume 
integral  (similar  to  equation  3.2)  at  a  radius  Rg .  This  always  leads  to  a  low 
frequency  spectral  peak  in  the  far  field.  Minster  (1973)  discussed  this 
problem  in  some  detail..  Snoke  (1976)  showed  that  this  approximation  leads 
to  acausal  results  in  the  time  domain  and  concluded  that  the  approximation 
was  incorrect  and  that  the  spectral  peaks  were  therefore  spurious.  Our 
results  are  somewhat  different  from  those  of  Archambeau 's  Rg  spectral 
peak.  That  approximation  removed  energy  from  the  low  frequency  part  of  the 
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spectrum  at  wavelengths  greater  than  R  .  Our  results  show  that  a  spectral 

s 

peak  can  in  fact  exist,  but  the  effect  is  to  add  in  additional  energy  at 
frequencies  slightly  above  the  "corner  frequency"  defined  by  the  pure  shear 
spectrum.  Also,  the  magnitude  of  the  spectral  peak  will  vary  with  angle. 

If  the  prestress  is  predominantly  pure  shear,  then  a  spectral  peak  will  be 
observed  only  near  the  quadrupole  nodes.  If  the  stress  were  more  inhomogeneous 
then  the  peak  can  be  observed  at  all  angles  of  observation.  Examples  of 
this  are  given  in  the  next  section.  There  have  been  a  number  of  observations 
of  spectral  peaks  from  earthquakes  and  explosions'  (e. g. ,  Linde  and  Sacks, 

1971)  and  an  inhomogeneous  prestress  is  a  likely  reason  for  this. 

Also  shown  in  figures  4,  5  and  6  are  the  waveforms  and  spectra  for  the 

relaxation/stress  pulse  term  (eq.  5.11)  alone.  Physically,  the  stress  pulse 
term  corresponds  to  a  sudden  stress  equal  in  magnitude  to  minus  the  initial 
prestress  on  a  spherical  surface  in  an  infinite  medium.  Thus  scattering 
diffraction  by  the  cavity  are  neglected  and  the  stress  pulse  is  rapidly 
radiated  away.  The  stress  pulse  term  is  a  fairly  good  approximation  to  the 
exact  solution,  especially  near  the  beginning  of  the  pulse.  The  main 
differences  are  that  the  stress  pulse  term  is  smaller  in  amplitude  and  stops 
abruptly  at  t  =  2R/v.  The  exact  waveform  looks  like  a  stretched  out  version 
of  the  stress  pulse  term.  The  spectra  contain  many  dips  not  present  in  the 
exact  solution.  These  are  due  to  the  abrupt  end  to  the  waveform.  These  are 
the  same  features  found  by  Burridge  (1975)  when  he  compared  his  exact  solution 
for  a  growing  spherical  cavity  in  a  pure  shear  field  against  Archambeau’s 
"transparent  source  approximation"  to  the  same  problem.  Use  of  the  stress 
pulse  term  alone  is  similar  to  the  transparent  source  approximation,  but  has 
seme  advantages.  The  transparent  source  required  the  solution  of  the  static 
problem  before  the  dynamic  problem  could  be  done.  Minster  (1973)  attempted  to 


use  the  transparent  source  approximation  for  an  ellipsoidal  rupture,  but  found 
the  static  problem  extremely  difficult.  Also,  Archambeau's  use  of  potentials 
have  some  singular  properties  which  lead  to  difficulties  in  other  problems 
than  the  case  of  pure  shear  he  considered.  In  particular,  for  the  case  of 
a  pure  compressive  stress  field,  the  method  gives  no  radiation.  The  stress 
pulse  term,  on  the  other  hand,  is  easily  evaluated  for  any  geometry.  It 
simply  requires  doing  the  surface  integral  (eq.  (3.3))  on  an  ellipsoidal 
surface  or  other  surface  of  interest.  It  is  singular  in  the  limiting  case 
of  the  dislocation  (a  flat  surface) ,  but  can  be  performed  for  any  surface 
enclosing  a  finite  volume. 
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6.  Preexisting  localized  stress  fields 

In  this  section  we  want  to  examine  two  physically  plausible  prestress 
fields.  The  first  is  due  to  a  point  source  of  compression,  the  second  due 
to  a  point  dislocation.  An  interesting  fact  is  that  both  of  these  stress 
concentrations  look  like  uniform  static  shear  field  when  far  away  from  the 
cavity,  but  have  some  anomalous  features  when  nearby. 

The  simplest  inhomogeneous  prestress  field  is  due  to  a  point  source  of 
compression.  Consider  a  spherical  region  in  the  earth  which  has  expanded  or 
contracted.  The  displacement  about  the  center  of  the  sphere  is  spherically 
symmetric  and  is  given  outside  the  spherical  region  by: 

u  =  -  Ve  =  cf;  (6.1) 

where  c  is  a  constant  depending  on  the  magnitude  of  the  contraction. 

We  now  consider  the  seismic  radiation  when  a  spherical  cavity  forms 
somewhere  outside  the  center  of  compression.  The  coordinates  are  shown  in 
figure  8.  To  translate  the  displacement  field  from  the  center  of  compression 
to  the  center  of  the  sphere  we  use  the  well  known  expansion  of  1/r  in 
spherical  harmonics. 

oo  £ 

7  -  £  V+T  Vcos  e)  R  <  L  (6.2) 

r  £=0  L* 

Taking  the  gradient  of  this  expression,  rewriting  in  terms  of  vector 
spherical  harmonics  and  using  equation  (5.8)  to  compute  the  normal  tractions, 


we  find 
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T(u  )  *e 


2yc  £ 


,£-2 


1-2  L 


£+1 


£(£-1)  P£0(e)  +  (£-1)  /£ (£+1)  B 


£0 


(6) 


(6.3) 


Note  that  for  L  >>  R 


T(u) 


A 

>e 


S3  -v  r 


(6. A) 


which  was  shown  in  section  4  to  be  a  pure  shear  stress  field,  as  noted  earlier. 

It  only  remains  to  compute  the  coefficients  of  ^u*.  The  actual  static 
displacement  field  cannot  be  expanded  as  in  the  last  section,  but  the  difference 
field  u*  can.  We  can  determine  the  coefficients  a.  ,  c.  by  comparison  of 
equation  (6.3)  with  equation  (5.9) 

For  L  »  R,  the  spectra  and  waveforms  are  identical  to  figures  2  and  3. 
Some  interesting  features  appear  when  the  cavity  is  closer  to  the  stress  con¬ 
centration.  Figure  9  shows  the  spectra  And  waveforms  near  the  quadrupole  nodes 
when  L  -  2R.  Instead  of  vanishing  at  these  angles,  the  amplitude  is  reduced 
by  about  a  factor  of  3.  The  spectra  of  both  the  P  and  S  waves  are  sharply 
peaked.  At  angles  away  from  the  nodes,  the  spectra  are  flat.  As  the  source 
is  moved  still  closer  to  the  cavity,  the  spectra  change  dramatically.  Figure 
10  show  the  spectra  near  the  quadrupole  maxima  when  the  stress  concentration 
is  .1  radius  from  the  cavity  surface.  The  comer  frequency  is  increased  by 
an  order  of  magnitude.  The  slope  above  the  comer  frequency  is  also  changed. 

The  corner  frequency  is  no  longer  a  measure  of  the  size  of  a  cavity.  Instead, 
it  is  a  measure  of  the  area  of  high  stress  on  the  cavity  surface.  An  important 
and  interesting  fact  is  that^'from  observations  of  first  motions  only,  the 
radiation  pattern  appears  to  be  very  nearly  that  of  a  quadrupole  source 
although  the  waveforms  and  spectra  vary  greatly  with  angle  (see  Figure  12). 
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As  a  second  example  we  consider  a  prestress  due  to  a  point  dislocation 
located  at  r  =  L.  The  coordinate  system  is  the  same  as  that  shown  in 
Figure  8.  The  initial  stress  field  and  the  spectra  and  waveforms  are  ^ 
considerably  more  complicated  than  for  the  point  source  of  compression,  but 
the  main  effects  are  the  same. 

The  displacement  field  due  to  a  dislocation  is  given  by  the  Volterra 
relation  (e.g.  Steketee,  1958) 


u(x)  =  /  u(x  ) 

•V  ^  /  */  •'vO 
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where  the  integral  is  over  the  surface  of  the  dislocation  and  G  is  the  static 
elastic  Green's  tensor.  Since  we  consider  only  a  point  Green's  tensor,  the 
integral  becomes: 


u (x)  =  S-T(G)-n 
x  ss 


(6.6) 


where  S  is  the  slip  vector  S  =  AuAA.  The  dot  products  with  n  and  S  are 
interchangeable.  The  static  Green's  tensor  is  given  by  Ben-Menahem  and  Singh 
(1968)  and  can  be  written  for  |x|  <  |x  | : 
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The  vectors  indexed  with  a  minus  are  the  same  as  those  in  equation  (5.1). 
The  remaining  eigenvectors  are: 


'  r*'1  t/TBETir  ♦  t  JtBI 


F+  _  rS.+l 
rvim 


i[(£+l+y)  /m+IT  /Blm  -  ( y-1-2 )  U+D^m]  (6*8) 
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If  we  restrict  the  problem  so  that  either  the  slip  vector  or  the  normal  n 

A 

is  in  the  direction  of  the  radial  vector  e^,  we  can  use  equations  (5.8)  for  the 
stress  derivatives.  We  take  the  slip  direction  to  be  radial,  take  S  5  |^sj,  and 
for  convenience  take  the  normal  vector  to  be  eg.  Then 


u(x)  -  S(T(G)-e  )-efl  (6.9) 

Most  of  the  terms  vanish  when  the  point  dislocation  is  located  at  6,4>  =0. 
Ben-Menahem  and  Singh  show  that  in  this  limit 
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After  some  tedious  but  straightforward  calculations,  we  find  that  the 
stress  field  at  the  surface  of  the  cavity  is: 
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where 


-  [nSt-i  +  s2'TS515' 


Sj^  -  A  (S.+1)  (2-y)  (2£+3)  +  2 p  (i+2-y)  (£+l)‘ 
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For  L  >>R  the  only  remaining  term  is 


T(u*)-ar  -  4j[(jtl  +  4-bu)  -  6  (-pt-i +  T-'t-J] 


(6. 


12) 


which  is  pure  shear  field  as  shown  in  section  4. 

Comparing  equation  (6.11)  with  equation  (5.9),  we  can  solve  for  the 

unknown  coefficients  an  ,  c.  ,  f.  .  We  then  use  equation  (5.11)  for 

Jom  Jem  Jem 

the  initial  value  field  and  use  the  technique  for  section  3  to  find  the 
exact  radiation  field. 

The  results  are  sh^$n  in  figure  11  for  a  dislocation  located  at 
L  =  1.5R,  one  half  radius  from  the  cavity  surface.  Both  the  P  and  S 
spectra  are  peaked  at  all  angles.  There  is  a  large  variation  in  the  slope 
above  the  "corner  frequency1.'  The  pulse  is  very  sharp  on  the  side  of  the 
sphere  closest  to  the  stress  concentration.  It  is  very  broad  on  the  other 
side.  The  dislocation  causes  a  region  of  high,  nonuniform  stress  on  the 
closest  part  of  the  cavity.  This  acts  as  the  source  of  the  radiation. 

The  pulses  are  widest  where  diffraction  is  most  obvious,  on  the  far  side 
of  the  cavity.  The  relaxation/stress  pulse  term  alone  matches  the  exact 
waveform  very  well  (except  for  the  diffracted  arrival)  for  angles  less 
than  90°,  but  rather  poorly  at  larger  angles.  This  is  to  be  expected 
since  the  surface  term  which  is  neglected  in  that  approximation  includes 
the  effects  of  diffraction.  An  interesting  feature  of  the  waveforms  is 
the  second  arrival  seen  at  low  angles.  This  is  a  diffracted  wave  which 
has  travelled  around  the  cavity  and  back  again.  An  analysis  of  the  first 
motions  shows  a  relatively  simple  angular  distribution.  The  distribution 
is  nearly  quadrupole  for  L  >  2R.  For  smaller  distances,  the  angular 
distribution  is  more  complex  (see  figure  12).  Note  that  the  first  motion 
for  angles  greater  than  90°  has  actually  reversed  direction.  The  dis¬ 
tribution  is  still  relatively  simple  and  could  easily  be  misinterpreted  as 
a  quadrupole  distribution  with  a  finite  number  of  observation  points. 
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Conclusion 

We  have  developed  a  general  method  for  computing  the  radiation  due  to 
the  instantaneous  creation  of  a  spherical  cavity  in  an  arbitrarily  prestressed 
medium.  Vector  spherical  harmonics  were  found  to  be  convenient  for  solving 
the  Green’s  tensor  integral  equation  in  the  frequency  domain.  Transforming 
the  initial  value  problem  to  a  boundary  value  (stress  pulse)  problem  proves 
to  be  a  great  simplification  and  provides  an  approximate  solution  as  well  as 
an  exact  solution.  The  approximate  solution  is  valid  when  diffraction  by  the 
source  is  not  too  important.  For  the  case  of  pure  shear,  the  results  of 
this  technique  agree  with  previous  solutions. 

We  have  examined  the  case  of  inhomogeneous  prestress  in  both  a  general 
manner  and  for  the  specific  cases  of  a  nearby  point  source  of  compression 
and  a  .dislocation.  Localized  stress  concentrations  result  in: 

1.  The  appearance  of  energy  above  the  usual  comer  frequency.  This 
can  increase  the  comer  frequency  substantially. 

2.  Far  field  spectra  peaks.  These  exist  near  the  quadrupole  nodes  for 
a  slightly  inhomogeneous  prestress  and  at  all  angles  for  a  very 
localized  stress  concentration. 

3.  Non-zero  amplitude  at  the  quadrupole  nodes. 

4.  The  creation  of  a  diffracted  wave  which  travels  around  the  cavity. 

5.  A  radiation  pattern  which  may  or  may  not  appear  to  be  quadrupole 
in  nature  but  is  generally  less  complex  than  the  change  in 
waveforms  and  spectra. 

6.  A  sharp  wave  on  the  side  of  the  cavity  nearest  the  stress  concen¬ 
tration  and  a  broad,  more  complex  wave  on  the  other  side. 
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The  results  given  here  may  be  useful  in  observational  studies  since  it  has 
implications  for  both  earthquakes  and  explosion  induced  tectonic  release. 

In  particular,  the  corner  frequency  is  often  used  as  a  measure  of  the  size  of 
the  source  R  s  v/u)c  where  R  is  the  cavity  radius  or  fault  length  depending 
on  the  type  of  source  being  studied  and  v  is  the  wave  velocity.  For  localized 
stress  concentrations,  the  corner  frequency  seems  to  give  an  estimate  of  the 
size  of  the  stress  concentration  which  may  be  much  smaller  than  the  source 
dimension.  Further,  the  existence  of  spectral  peaks  could  lead  to  errors  in 
estimates  of  seismic  moments.  The  seismic  moment  is  defined  by  M  ■  yuA  where 
y  is  the  shear  modulus,  u  is  the  mean  displacement  and  A  is  the  surface 
area  of  the  ^lilt.  It  is  also  proportional  to  the  zero  frequency  far  field 
spectral  amplitude  (e.g.,  Kostrov  (1974)).  If  observations  are  made  of 
spectra  only  close  to  the  corner  frequency,  the  apparent  low  frequency  amplitude 
could  be  much  higher  than  the  actual  zero  frequency  amplitude  leading  to 
overestimates  of  seismic  moment.  This  could  lead  to  misinterpretation 
involving  not  only  fault  length,  but  stress  drops  as  well.  In  addition  to 
these  problems,  great  care  must  be  taken  to  remove  the  effects  of  the  source 
when  doing  studies  of  earth  structure  since  the  variation  of  pulse  shape 
with  ("take-off")  angle  could  cause  confusion  between  the  effect  of  earth 
structure  and  the  effect  of  localized  stress  fields. 

On  the  other  hand,  these  results  should  make  it  possible  to  study  the 
local  stress  variations  in  the  earth.  The  techniques  used  in  section  5  can 
be  used  to  invert  for  the  stress  (difference)  field  after  an  explosion.  If  a 
large  earthquake  is  to  occur  in  some  region,  it  means  that  a  large  stress 
concentration  has  developed  which  is  not  relieved  by  small  earthquakes  nearby. 
Observations  of  a  change  in  the  waveforms  and  spectra  of  small  earthquakes  in 
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an  area  could  indicate  a  developing  stress  concentration.  Such  an  observation 

V/J 

was  made  recently  by  Zolljjeg  (1979),  who  reported  a  change  in  the  high  frequency 

slope  of  the  spectra  of  small  earthquakes  in  Missouri  from  w  ^  before  a 

_2 

large  earthquake  to  w  after  the  earthquake.  He  also  reported  a  variation 
of  the  slope  with  ^STj-muth.  This  is  the  sort  of  behavior  that  is  to  be  expected 
if  a  local  stress  concentration  is  relieved  by  the  large  earthquake.  In 
addition  to  this,  observations  of  the  features  described  in  this  paper  in 
the  variation  of  spectra  and  waveforms  with  angle  could  be  used  to  identify  the 
location  of  developing  stress  concentrations. 
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Figure  Captions . 


Fig.  1. 


Fig.  2. 


Fig.  3. 


Fig.  4. 


Fig.  5. 


Fig.  6. 


The  problem  considered  here.  A  cavity  suddenly  forms  in  a  prestressed 
medium.  The  radiation  due  to  the  relaxation  of  the  stressed  material 
can  be  computed  using  a  Green's  tensor  integral  equation  with  the 
aid  of  vector  harmonics  appropriate  for  the  cavity  geometry. 

Far  field  P  and  S  wave  spectra  for  a  uniform  shear  field.  The  first 
figures  show  the  spectra  for  the  stress  pulse  term  alone.  The 
second  are  the  exact  solutions.  The  dips  in  the  spectra  are  removed. 
The  frequency  is  in  units  of  toR/V.  All  examples  in  this  paper  use 
o=8  km/sec,  g  =  5  km/sec. 

Far  field  pulse  due  to  the  creation  of  a  stress  free  sphere  in  a 
uniform  shear  field.  The  exact  solutions  agree  with  the  solutions 
of  Hirasawa  and  Sato. 

The  most  general  Z  *  2  mode  radiation  consists  of  three  independent 
waves-two  coupled  P  and  S  waves  labeled  PI,  SV1  and  P2,  SV2  and  an 
uncoupled  S  wave  labeled  SH.  The  spectra  for  the  Z  =  2  modes  are 
shown  here.  Two  modes  are  flat  and  three  are  peaked.  All  higher 
modes  (£  >  2)  have  peaked  spectra. 

Waveforms  for  general  Z  =  2  case.  A  particular  linear  combination 
of  these  produces  the  pure  shear  waveforms  shown  in  Figure  3.  The 
second  column  shows  waveforms  using  the  stress  pulse  term  alone.  The 
exact  solution  is  a  stretched  out  version  of  the  approximate  solution. 
Waveform  for.  Z  =  4.  All  higher  l  modes  are  oscillatory.  The 
number  of  oscillations  increases  with  £. 


51 


Fig.  7. 


Fig.  8. 

Fig.  9. 


Fig.  10. 


Fig.  11. 


The  effect  of  an  inhomogeneous  prestress  is  to  add  in  energy  at 
frequencies  higher  than  the  usual  corner  frequency  corresponding 
to  uniform  shear.  Shown  here  are  the  spectra  of  the  far-field 
first  P-wave  for  £  =  3  to  £  =  6  and  the  pure  shear  spectrum  for 
comparison.  A  sufficiently  inhomogeneous  prestress  can  result  in 
a  low  frequency  spectral  peak  which  varies  in  magnitude  with 
angle. 

Coordinate  system  used  when  a  cavity  is  created  near  a  preexisting 
center  of  compression. 

Far  field  radiation  near  the  quadrupole  nodes.  When  the  center  of 
compression  is  located  one  radins  away  from  the  cavity,  the  radiation 
field  is  much  like  a  pure  shear  field.  Near  the  quadrupole  nodes 
however,  there  is  a  substantial  difference.  The  displacement  does 
not  vanish.  It  is  reduced  from  the  maximum  by  about  a  factor  of  3. 
The  pulse  is  oscillatory  and  the  spectrum  is  peaked. 

Far  field  spectra  near  quadrupole  maxima.  The  spectrum  is  strongly 
affected  by  the  location  of  the  stress  concentration.  Note  the 
increase  in  corner  frequency  and  the  change  in  slope  near  the  corner. 
Radiation  field  from  spherical  cavity  created  one-half  radius  from 
a  point  dislocation.  Shown  here  are  the  far  field  spectra  and 
waveforms  as  a  function  of  angle  around  the  cavity.  The  spectra  are 
all  strongly  peaked.  The  waveforms  are  narrow  on  the  side  of  the 
sphere  near  the  stress  concentration,  broad  on  the  other  side.  The 
second  arrival  seen  on  the  first  three  shear  waves  is  a  diffracted 


phase  which  has  travelled  around  the  cavity. 
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Fig.  12.  First  motions  for  localized  prestress  fields.  For  L  >  2R  ,  the 
first  motions  are  very  nearly  quadrupole  in  nature.  At  closer 
distances,  the  patterns  shift.  For  the  center  of  compression,  there 
is  simply  a  change  in  the  "nodal"  angles.  For  the  dislocation 
a  pulse  directed  opposite  to  the  quadrupole  pulse  gradually  becomes 
prominent  at  large  angles.  The  angular  distribution  of  first  motions 
is  less  complex  than  the  angular  variation  of  spectra  and  waveforms. 
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The  Synthesis  of  Complete  Seismograms 
In  an  Earth  Model  Specified  by  Radially 
Inhomogeneous  Layers 


Vernon  F.  Cormier 
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ABSTRACT 


Traditional  methods  of  synthesizing  seismograms  in  the  near 
or  regional  distance  range  (0  to  500  km.)  describe  the  earth  model 
with  planar  homogeneous  layers.  In  the  high  frequency  band  (.2  - 
10  Hz.),  in  which  the  seismic  ground  motion  is  primarily  observed, 
uniformly  asymptotic  solutions  to  the  depth  eigenfunctions  can 
instead  allow  a  radially  symmetric  earth  model  to  be  described 
by  inhomogeneous  spherical  layers.  Ground  displacement  u  is  cal¬ 
culated  in  the  frequency  domain  by 


r 


where  r  is  a  contour  in  the  complex  ray  parameter  (p)  plane,  and 
M(u>,p,  4>,  A)  a  point  representation  of  the  earthquake  or  explosion 
source  including  the  effect  of  horizontal  propagation  to  A. 

The  response  of  the  earth  model  f(u,p),  calculated  from  the 
propagator  matrix  equation  for  a  source  in  a  radially  inhomo¬ 
geneous  sphere,  includes  all  possible  body  and  surface  waves.  Airy 
functions  are  chosen  to  define  the  inhomogeneous  layer  matrices. 
Numerical  difficulties  usually  encountered  in  the  calculation 
of  f(ui,p)  in  layered  media  are  avoided  by  the  calculation  of  sub¬ 
determinants  of  the  fundamental  matrix  solution  and  by  the  de¬ 
composition  of  the  propagator  matrix  in  a  layer  into  a  sum 
of  matrices  of  differing  numerical  order  whenever  the  Airy  functions 
behave  exponentially.  The  contour  integral  can  be  evaluated 
by  the  residue  theorem  or  by  numerical  integration,  the  time  domain 
response  obtained  by  fast  Fourier  transform.  Most  applications 
reguire  an  earth  model  to  be  described  by  no  more  than  four  to 
five  inhomogeneous  layers. 
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INTRODUCTION 


At  teleseismic  distances  in  the  frequency  bend  (0.2  -  10  Hz.), 
surface  and  direct  body  waves  are  sufficiently  separated  in  time 
that  either  may  independently  be  used  to  extract  source  and  struc¬ 
tural  properties.  In  the  distance  range  0  to  500  km.,  however, 
a  complete  or  much  larger  set  of  seismic  waves  must  be  employed 
to  properly  predict  ground  motion  in  even  a  small  time  window. 
Although  complicating  the  problem  of  waveform  synthesis,  the  short¬ 
er  signal  duration  offers  an  opportunity  with  high  frequency 
recordings  for  resolving  finer  structural  details  and  near  field 
source  properties  impossible  with  longer  period  teleseismic  re¬ 
cordings  . 

At  sufficiently  high  frequencies  asymptotic  solutions  for 
the  radial  or  depth  eigenfunctions  accurately  approximate  the  exact 
solutions  for  radially  inhomogeneous  layers  even  with  relatively 
severe  velocity  gradients.  Formalisms  developed  by  Langer  (1932, 
1949)  ,  Olver  (1954)  ,  and  VJasow  (1965)  can  be  employed  to  obtain 
asymptotic  solutions  (in  frequency)  to  either  the  radial  component 
of  seismic  potentials  (Richards,  1976)  or  to  the  radial  vector- 
matrix  equations  of  seismic  motion  (Chapman,  1975;  Woodhouse, 

1978) .  Such  solutions  allow  a  much  simpler  description  of  com¬ 
plicated  earth  models  by  a  small  number  of  radially  inhomogeneous 
layers.  By  their  validity  at  high  frequency  these  solutions  are 
particularly  suited  to  applications  in  the  frequency  band  in  which 
near  and  regional  seismic  ground  motion  is  primarily  observed  and 
interpreted.  Compared  to  a  model  specified  by  planar  homogeneous 
layers,  a  model  specified  by  radially  inhomogeneous  layers  will 
usually  require  many  fewer  mathematical  operations  to  construct 
its  frequency-ray  parameter  response.  This  objective  is  met 
in  the  following  sections  by  detailing  a  procedure  for  synthe¬ 
sizing  complete  seismograms  in  such  an  earth  model.  The  pro¬ 
cedure  combines  (1)  the  zeroth  order  (in  frequency)  asymptotic 
solutions  to  the  propagator  and  fundamental  matrix,  (2)  the 
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notation  of  generalized  vertical  slownesses  or  cosines  de¬ 
veloped  by  Richards  (1976) ,  and  (3)  the  vector-matrix  methods 
of  Abo-Zena  (1979)  for  avoiding  the  numerical  difficulties  in 
calculating  the  response  of  a  layered  medium.  The  THEORY 
section  of  the  paper  develops  notation  and  the  generalized 
solution  for  displacement.  The  reader  primarily  interested  in 
developing  practical  numerical  codes  may  proceed  directly  to  the 
section  titled  EXAMPLE  USES  OF  GENERAL  SOLUTION  FORM,  referring 
as  needed  to  the  sections  EVALUATION  OF  PROPAGATOR  PRODUCTS , 
THEORY',  and  Appendices  A  and  B. 

THEORY 

General  Solution  Form 

By  applying  the  vector  representation  theorem  for  a  sphere, 
the  seismic  displacement  vector  u  =  u^r  +  u„S_  =  may  be 
written  as 

/to  -i.  tiit  °°  n  m  m  m 

due  I  1  (UP  +  VB  +  WC  )  (1) 

-«>  n=0  m=-n  —n  ~ n  ~ n 

iti  rci  in 

where  £  ,  B  ,  and  C  are  vector  spherical  harmonics  defined  by 
n  n  n 

m  m 

P  =  rY  (0,$) 

n  n 

B1"  =  V1Ym(0,4;)/[n(n  +  1)]  **  (2a) 

~ n  n  L  J 

C  =  -rxVjY  (  0  ,  4> )  /  [n  ( n  +  1)]  * 
n  n 

=  i  Jq  +  coseceo  |^)  . 

mm  m 

Here  the  P  ,  B  ,  and  C 
n  n  ~ n 


are  fully  normalized: 
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m 

Y  (6,$) 


n 


(-) 


m 


j"2n+l  (n-m) 


L  4  71  (n-m) 


t] 


m 

P  (cos  6) 
n 


iirnp 


sin6ded<i 


sin6d6dg 


sinSdBco 


=  6mm'  6nn '  . 


(2b) 


(2) 


£,  6_,  and  $  are  unit  vectors  in  the  spherical  coordinate 

m 

airections  r,6,<J>.  P  (cos6)  is  an  associated  Legendre  function. 

n 

The  radial  eigenfunctions,  U,  V,  and  W  satisfy  vector-matrix 
equations  of  the  form 

VT  =  O)(A0  +  U)  1 A  x  T  +  Cl)  2A2^)v^  (4) 


for  SK  type  body  waves  and  toroidal  free  oscillations,  and 

-2*  S.  TS  . 


VS  =  Cl)  (A  S  +  U)-3A,S  +  Ci)~2  A  S>)vt 
=0  =1  —2  — 


8r  - 


5) 


for  P-SV  type  body  waves  and  spheroidal  free  oscillations  (e.g., 

Chapman,  1973).  When  gravitational  effects  can  be  ignored, 

T  S 

the  displacement-stress  vectors  v  ,  v  are  defined 


(6) 


where  the  radial  stress  functions,  T,  R,  and  S  are  defined  in 

terms  of  U,  V,  and  W  by  Alterman  et  al.  (1959) .  Here  and  later, 

quantities  written  vertically  within  brackets  denote  a  column 

vector  and  quantities  written  horizontally  and  separated  by 

commas  denote  a  row  vector. 

A  matrix  F  satisfying  either  (4)  or  (5)  in  the  column 
T  s” 

vectors  v  ,  v  is  defined  as  a  fundamental  matrix  (Gilbert  and 
Backus,  1966)  or  matricant  (Gantmacher,  19T5)  .  As  shown  by 
Gilbert  and  Backus  (1966),  the  fundamental  matrices  of  these 
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equations  allow  a  convenient  and  notationally  co-pact  method 
of  satisfying  the  boundary  conditions  at  layers  cf  an  earth 
model  having  discontinuities  in  elastic  moduli  and/or  density. 
The  stress  displacement  vector  v  at  a  radius  "a"  is  related  to 
that  at  a  radius  rn  through  a  product  of  propagator  matrices 
constructed  from  the  fundamental  matrices  for  each  i-th 
layer : 


v(a)  =  K (a,r  ) v (r  )  =  K  ...  K  v (r  )  , 
—  —  n  n  *=  i  — n —  n 


where  K^r.^,  r±)  =  £i(ri-l)£i~  {ri} 


with  r^_^  and  r^  denoting  the  radii  bounding  the  top  and 
bottom  boundaries  of  the  i-th  layer  respectively.  Figure  1 
illustrates  the  layering  scheme  and  conventions  used  for  indices. 

The  general  solution  to  (4)  or  (5)  for  the  displacement- 
stress  vector  v,  including  the  effects  of  a  source  vector 
term  y,  can  be  written  as 


v  (a)  *  S  (a,r0 )  v(ro )  +  K(a,C)v(Uc;  (9) 

o 

where  ro  is  some  reference  radius  at  which  a  starting  solution 
v(ro)  is  defined  (Wasow,  1965).  The  complete  displacement  solu¬ 
tion  can  then  be  obtained  by  substituting  the  results  for  U,V, 
and  W  determined  from  (9)  into  the  displacement  representation 
given  by  (1) . 


Fundamental  Matrices  of  Radially  Inhomogeneous  Lavers 


The  basic  form  of  solution  for  the  disolacement-stress 


vector  given  by  (9)  is  independent  of  any  representation  of  the 
earth  model,  e.g.,  planar  homogeneous  layers,  spherical  inhomo¬ 
geneous  layers,  etc.  Assume  now  that  the  earth  model  is  speci¬ 


fied  by  radially  inhomogeneous  layers  and  apply  the  zeroth  order 
uniformly  asymptotic  solutions  to  the  fundamental  matrix  and 
propagator  of  such  layers.  Adopting  the  notation  and  solutions 
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given  by  Aki  and  Richards  (1980),  omitting  scarce  normalization , 
and  rearranging  rows  and  columns  to  agree  with  the  stress- 
displacement  vector  used  by  Woodhouse  (1978)  gives 


r(r)  =  r-Sf 


F  (r) 


-  i 


\P 


C5(1) 

-ir-1pg(1) 

-Vg(3) 

-ir-»g(3) 

-  ■  i  n 

r” 1 ph  ^ 3  ^ ■ 

-iYh(3) 

-xAg(1) 

-xAg(3) 

iBnh(1) 

— ;snh(3) 

BCg(1) 

RV  (3) 

-B^g 

Ah(1' 

Ah(3) 

(10a) 

xAg<3) 

BYg(3) 

-Yg:3) 

.  - 1 

-  x  r  pc 

-xAg{1> 

B£g{1) 

-£S(1) 

ir~ " pg ^ 

-xBnh(3) 

Ah(3> 

-r-:?h(3) 

*Th(3> 

--xBnh(1) 

-Ah(1) 

r_1nh(1) 

xYh(1> 

(3) 

(3) 


(10b) 


for  P-SV  waves  where  A=2r-2p2 y-p ,  B=2r~  1  p  ,  x=’/-T,  y  is  the  shear 
modulus  profile  y(r)  and  p  the  density  profile  p  )r) ;  and 


F*(r)  =  r_1/x 


Frl(r)  =  r^f 


y_l5h(1) 

iyJsnh(1) 

-xy%<3> 

-iy^h(1) 


y"V3>  ' 

-iy%(3) 

,-vf3) 

-y  h 

-%•_  (D 
y  ^ 
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for  SH  waves,  g^,  g^3\  h^,  h^  are  generalized  radial 

eigenfunctions  defined  in  terms  of  Airy  functions  (see  Appendix)  . 
,  v  ,  v 

£,  £r  hf  n  are  generalized  vertical  slownesses  defined  with 
the  generalized  radial  eigenfunctions  by 

.  „.(D 


£  - 


,  (1) 


u>xh 


(1) 


v  _a.(3) 
£  =  -2- 


uig13' 


(12) 


Y.  zH 


,  (3) 


uih 


(2) 


M* 


The  ray  parameter  p  in  the  fundamental  matrices  can  be 
associated  with  the  degree  n  of  the  vector  spherical  harmonics 
in  equation  (1)  by 

p  =  (n  +  h)/u.  (13) 

Functions  analytic  in  radius  specify  the  seismic  velocity 
behavior  in  all  layers.  Appendix  A  outlines  and  references 
methods  of  fast  evaluation  of  the  radial  eigenfunctions  and 
generalized  vertical  slownesses  in  such  profiles. 

The  fundamental  matrices  may  be  defined  using  any  one  of 

several  different  pairs  of  linearly  independent  Airy  functions. 

The  pair  Ai(-z)  and  Bi(-z)  chosen  by  Woodhouse  (1978)  represent 

standing  waves  that  exponentially  decay  or  grow  with  increasing 

i2v/3 

depth  below  a  ray  turning  point.  The  pair  Ai(-ze  )  and 

Ai(-z  "  ''  )  chosen  by  Richards  (1976)  represent  up-  or  down¬ 

going  travelling  waves.  The  intermediate  choice  of  Ai  (-ze1'271^3) 
and  Ai(-z)  still  allows  identification  of  travelling  waves 
within  a  single  layer,  but,  as  will  be  shown  in  a  later  section, 
avoids  a  numerical  problem  in  the  computation  of  the  response 
of  a  layered  model.  The  superscript  d)  in  equations  (10a-12) 
refers  to  the  use  of  Airy  function  appropriate  for  upgoing 
travelling  waves,  Ai  (-z  ~2Ti//3)  .  The  individual  layer  propagator 
defined  using  the  set  Ai  (-ze  *'^T7/^)  ,  Ai  (-ze  equals  that 

defined  using  the  set  Ai  (-ze^2  ,  Ai(-z).  Thus,  provided 

the  propagator  wuthin  a  layer  is  formed  by  (8),  the  superscript 
(3)  can  refer  to  either  the  use  of  Ai  (-ze-1'2  "^J)  or  the  use 
of  Ai(-z)  when  convenient. 

Boundary  Conditions 

The  propagator  formalism  satisfies  boundary  conditions  on 
stress  and  displacement  components  at  each  layer  interface.  The 
only  remaining  boundary  conditions  are  the  vanishing  of  stress 
at  the  free  surface  and  a  radiation  condition  in  the  lowermost 
layer.  For  a  sphere  of  radius  "a"  the  free  surface  condition 
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requires  that  R(a) =£ (a) =T (a) =0 .  Two  possible  boundary  conditions 
may  be  invoked  in  the  last  layer.  (1)  The  last  layer  possesses 
only  standing  waves  that  exponentially  decay  with  increasing 
depth  below  their  turning  points.  (2)  If  it  is  desired  to 
exclude  any  waves  from  propagating  in  the  last  layer  from  the 
solution  at  the  free  surface,  the  last  layer  must  possess  only 
downward  propagating  waves.  Henceforth  the  term  "decay" 
will  refer  to  boundary  condition  (1)  and  "reflectivity"  to 
boundary  condition  (2) .  Witn  these  boundary  conditions  and 
when  the  superscript  (3)  refers  to  the  Airy  function  Ai(-z), 
the  stress-displacement  vector  for  P-SV  waves  has  the  form 

for  decay 


for  reflectivity. 


or 


F(r„) 


„) 


0 
b 
0 
b. 

b. 


For  SK  waves  the  form  becomes 


v(r0)  =  F'(r0) 
or  =  F‘(r  0 ) 


for  decay 


for  reflectivity. 


The  quantities  b2,  b2 ,  and  b3  are  constants. 

The  analytic  velocity  profiles  in  the  last  inhomogeneous 
layer  specify  the  velocity  up  to  the  center  of  the  earth. 

They  need  only  accurately  model  the  P  and  S  velocities,  however, 
up  to  the  maximum  depth  that  rays  bottom  in  the  ray  parameter  (p) 
domain  needed  for  the  solution  of  a  particular  problem.  The 
velocity  profiles  may  have  multiple  turning  points  but  only 
one  turning  point  may  exist  in  this  p  domain.  Because  radial 


75 


inhomogeneitv  and  sphericity  allow  rays  to  have  turning  points, 
the  last  layer  automatically  radiates  seismic  energy  back 
towards  the  surface  with  the  decay  boundary  condition.  With 
the  reflectivity  condition  applied  no  the  last  layer,  only 
waves  reflecting  from  or  bottoming  above  the  last  layer  reach 
the  surface. 

The  reference  radius  ro  can  be  imagined  to  be  far  below 
the  lowest  level  for  which  rays  turn  in  the  p  domain  of  interest. 
A  particular  rp  need  not  be  specified  because  the  functions 
evaluated  in  r0  cancel  in  the  subsequent  solution  procedure. 

Including  the  free  surface  and  decay  boundary  conditions 
in  the  general  solution  form  of  (9)  gives 


U(a)~ 

'0  ' 

V(a) 

=  K  ( a ,  r  )F  (r  .) 

b, 

0 

=  n  =n  n-1 

0 

0  . 

_b2_ 

f&  K(a,UY  U)d? 

J  r. 


(14) 


for  P-SV  waves,  and 


W  (a) 
0 


-  £'(a'r„>S'nlrn-l>  [£,]+/a  S'(a.5)v-(C)« 


(15) 


for  SH  waves.  The  unprimed  and  primed  fundamental  matrices, 
propagators,  and  source  vectors  correspond  to  the  ?-£7  and  SH 
cases  respectively. 

Solution  in  a  Layered  Sphere 

Multiplying  (14)  by  (r  ^) K (rn , a)  and  (15)  by 

F'  1 (r  ,)K'{r  .a)  and  rearrancing  terms  cives 
=  n  n-i  —  n 
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Lb2J 
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for  P-SV  waves,  and 
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for  SH  waves,  where 
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=  (rn-l,a)  _  ^n-l^n-l^n^  *  *  *  =  i  (r  i ,a) 
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(20) 


R’  J,  S’  1  are  defined  analogously  from  the  primed  fundamental 
matrices  and  propagators.  The  sbuscript  s  refers  to  the  layer 
index  of  the  source  layer,  rg  denoting  the  radius  at  its  lower 
boundary  (see  Figure  1) . 

Solutions  of  equations  (16-17)  for  the  radial  functions 
U(a),  V(a)  and  W(a)‘  are 
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V  (a)  = 


- 1 

G  R 
3  3  1 


-  1  -  1 
R  1  3  R  3  2 


G3R3! 
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(22) 


W  (a)  = 


G' 


R! 


-  j 


i  i 


(23) 


-  l  - 1 

where  R. .  denotes  the  rj-th  element  of  R 
3-1  J  = 

EVALUATION  OF  PROPAGATOR  PRODUCTS 

Care  must  be  taken  in  the  evaluation  of  the  propagator 
products  of  equation  (20)  not  to  lose  precision  from  the  sub¬ 
traction  of  large  quantities  of  nearly  equal  value.  Knopoff 
(1964),  Dunkin  (1965),  and  recently  Abo-Zena  (1979)  have  shown 
how  to  avoid  this  difficulty  in  models  having  planar  homogeneous 
layers.  All  of  these  methods  propagate  a  matrix  of  minors  of 
the  inverse  fundamental  matrix.  The  methods  used  by  Abo-Zena, 
however,  free  the  layered  model  from  any  constraints  on  its 
total  thickness.  Because  the  precision  problem  also  exists 
for  models  composed  of  radially  inhomogeneous  layers,  the  matrix 
methods  developed  by  these  authors  can  similarly  be  applied. 


The  R. .  Elements 
3-3 _ 

The  elements  R. 1  and  R.1  in  (21-23)  may  be  written  as 
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where  the  row  vectors  Ep,  EgV,  EgH  are  deterzr.ired  by 
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Because  the  elements  r7^  do  not  depend  on  tie  r 

of  the  spherical  harmonics,  they  may  be  moved  c: 

m  summation  when  (21-23)  are  substituted  in  equa 

is  then  convenient  to  perform  the  summation  ever 
m  m  m  m  m 

products  G,P  ,  G  P  ,  G  B  ,  GB  .  and  G'C  .  The- 
^  1  — n  3  — n  i  — n  3  — n  i  — n 

may  be  written  as: 


der  number  m 
side  of  the 
ion  (1) .  It 
m  of  the 
quantities 
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Appendix  B  defines  the  column  vectors  s_r ,  s_- ,  s__ ,  s_' ^  and 
for  generalized  earthquake  and  explosion  point  sources, 
and  t'  are  the  sums  of  the  column  vectors  defined  in  equation 
(B7 )  in  the  vector  directions  r,0,  and  6  : 


1  =  £rr  +  se6  +  s^i 

i’  =  s’ee  +  s’^i 


(27) 


The  Displacement  Representation 

The  theorem  for  the  transpose  of  the  product  of  matrices  AE 
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and  the  fact  that  the  transpose  of  a  scalar  leaves  it  unchanged, 
imply  the  identities 
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The  transposes  of  t  and  t'  are  taken  as 
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£r  etc.  are  now  row  vectors.  Using  equations  (1) ,  (21-24) 
and  (29)  ,  it  then  follows  that  for  a  generalized  point  source 
and  for  a  model  described  by  radially  inhomogeneous  layers, 
the  displacement  u  can  be  written  as 


u(r,8,$,t)  = 

OO  CD 


tTKT[EpTEsv  -  t*,  Ep)K 

i  i 

O  O  (->  o 

l _ 1 

[l,0,0f0]KT[EpT  Esv  -  EsTvEr 

,3K 

0 

1 

0 

-  .  i 

0 

— SH— ' 


wdpdio 


(31) 


The  n  summation  has  been  converted  to  an  integral  over  p  by 
taking  the  first  term  in  a  Poisson  transform  (e.g. ,  Nussenzweig, 
1965)  . 


Matrix  Multiplication  Scheme 

T  T 

The  matrix  [Ep  Egv  -  Egv  E  ]  is  an  anti-symmetric 
matrix  of  minors  of  F  “ 1  of  the  form 


(cont 'd) 


(32) 
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‘n-1 


n-1 


„(3).  (3) 
g  h 


where 


d2  =  -(A2  +  b2^K) 

d2  =  (r  *pA  +  Bn£) 

,  •  v 

d3  =  tp n 

d„  =  ipC 

s/  y 

d5  =  (n£  +  r  2p2) 


(32) 

(cont 


(33) 


If  the  reflectivity  instead  of  the  decay  boundary  condition  had 
been  used  in  the  last  layer,  the  Airy  function  Ai  (-ze-^^11^3) 
would  replace  Ai(-z)  in  the  definitions  of  the  vertical  slow¬ 
nesses  appearing  in  the  quantities  d^. 

All  elements  of  Yfi  are  evaluated  at  the  boundary  of  the 
last  layer  and  in  the  analytic  velocity  profiles  of  the  last 
layer.  Evaluation  of  the  P-SV  terms  in  the  brackets  of  (31) 
involves  a  successive  redefinition  of  a  matrix  at  each 
i-th  higher  layer  boundary  starting  at  the  lowest  layer  boundary: 


-n-1 


Sn-lMn-1 


Y 

*=n- 


2 


T 

K  -Y  ,  K 
— n-2-=n-l— n- 


2 


Xi 


I2L 


1 


Each  newly  defined  Y^  matrix  remains  an  antisymmetric  matrix. 
Only  five  independent  terms  are  needed  to  specify  each  new 
(Menke,  personal  communication). 
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The  denominator  for  the  P-SV  term  is  given  by  the  inter¬ 
section  of  the  first  row  and  the  second  column  of  the  result 
for  Yi .  The  numerator  is  calculated  by  first  evaluating  the 
product 


v 

“S~ 


1 


KYnK 


(35) 


right  multiplying  by  Kg (r  ,rg_^) - . .K x  (r j ,a) ,  and  left 

multiplying  the  second  column  of  this  result  by  the  row 
T 

vectors  of  t  . 

The  methods  of  Abo-Zena  (1979)  can  be  applied  to  evaluate 
the  products  of  equation  (31)  for  the  P-SV  term  and  the  products 
E-hK ' ,  E_„K '  for  the  SH  term.  These  methods  decompose  the 
propagator  matrix  of  a  layer  into  a  sum  of  matrices  of  differing 
numerical  order  whenever  the  depth  eigenfunctions  are  exponentially 
small  or  large.  In  an  inhomogeneous  layer  this  happens  at  ray 
parameters  for  which  a  P  or  S  wave  bottoms  far  above  the  layer. 

At  such  ray  parameters  only  tunneled  energy  can  reside  in  the 
layer.  Because  tunneling  is  a  frequency  dependent  phenomenon, 
a  particular  exponential  decay  value  depends  on  frequency. 

When  this  situation  occurs  in  the  i-th  layer,  define  the 
positive  integers  iru  ( j  =  1  —  4 )  and  the  complex  numbers  (j  =  l-4) 
with  magnitude  of  order  1  from  the  products 


g{1)  (ri)g(3)  (r.^) 

= 

^emi  , 

g(3>(ri)g^Nr..1) 

= 

S2e‘m2, 

(36) 

h(1>  (ri)h(3) (ri_1) 

= 

4,6®*  , 

n(3)  (ri)h(1)  (ri_1) 

* 

The  Sj  and  m^  values  carl  be  calculated  and  returned  by  an  Airy 
function  subroutine  based  on  the  exact  and  asymptotic  formulae 
given  in  Abramowitz  and  Stegun  (1964).  (More  detailed  references 
are  given  in  Appendix  A.)  Next  define  the  matrices  Vjh j  from 
the  product  of  column  vectors  v ^  and  row  vectors  h ^ .  After  the 
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exponential  scale  factors  of  the  fundamental  matrices  are 

removed,  v ^  is  given  by  the  jth  column  cf  the  F_.  (r^)  matrix  and 

h.  by  the  jth  row  cf  the  F. (r.  ,)  matrix.  The  crooacator  K. 

— j  J  =ii-l  -  '  =i 

can  then  be  written  as 


K .  = 
=1 


l2Ji2  e 


-m2 


X3^3e‘ 


-n  l 


(37) 


In  Abo-Zena  (1979)  such  a  decomposition,  of  is  invoked 
to  varying  degrees  depending  on  which  of  three  diferent 
relations  exist  between  phase  velocity  and  layer  velocities. 
For  inhomogeneous  layers  these  cases  are  restated  in  terms  of 
the  exponential  behavior  of  the  Airy  functions  used  in  the 
definition  of  the  fundamental  matrices: 


I.  m,  =  n,  =  m,  =  m„  =  0 

1  i  3  »♦ 


All  radial  eigenfunctions  behave  either  as  sinusoidal  or  phasor 
type  functions  having  a  magnitude  of  order  I.  Compute  directly 


and  determine  Y.  from 


-  £iTIi*i£i 


II.  m,  and  m  =  0,  m,  and/or  m. ■ 

3  **  i  A 


A  matrix  oc;  may  be  defined  as 


22  “  hki  +  v4h, 


(36) 


Compute  Yi  from 

Si  =  <22>TXi+122  +.  ^:TX!T+  *2  2? 

+  err'l|[h!T  Vj1  +  (o£)  T]  Yi+1  (Vjhj+  oc)  -  (££)TYi+1££j  (39) 

5_mj|[h2Tv2T  +  (££) T]  Vi+1  (v,h2  +  oc)  -  (oc)  TYiflog j 
III.  m  and/or  m2  ^  0 ,  m3  and/or  mt  i-  0 


+  e 
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Compute  Y.  from 

rn»T»n,'rp  rprprprp 

li  =  +  *2i2>2i+i^2+  z,!*,)  +  <*>3+  h;v;)Yi+1(Lht  +  v3h3) 

+  .«»**•>  .a,  *  v,h.) 

+  e  (mi  -ms )  (^T  +  hVly.^lv^  +  v.h,)  (40) 

+  e  (ms-mi )  ,,££  +  h^)Y.+1(v,h.  +  v2h2) 

+  e"(mi+in3)  (hV  +  hV)Y.  ,  (v  h  +  v  h  ) 

—2—2  —4—1. '=1  +  1  —4—1.  —2—2 


All  of  the  symmetries  noted  by  Abo-Zena  can  be  exploited  in 
the  evaluation  of  the  matrix  products  of  algorithms  I-III. 

The  propagator  products  in  the  SH  calculation  can  be 
evaluated  using  simpler  analogous  algorithms  for  the  decomp¬ 
osition  of  the  SH  propagator.  New  row  vectors  ESH  are  given 
by  right  multiplying  by  the  propagators. 

At  large  iru  the  product  having  the  largest  exponential 
factor  in  algorithms  II  and  III  dominates  the  final  result 
for  and  the  other  terms  may  be  ignored.  The  exponential 
behavior  assumed  here  is  based  on  the  travelling-standing  choice 
for  the  Airy  functions  used  in  defining  the  fundamental  matrices, 
Ai(-ze  *'^7T//3)  and  Ai(-z).  The  displacement  representation 
given  by  (31)  equals  that  given  by  substituting  Ai(-ze  ^27T/^) 
for  Ai(-z)  in  (31)  in  all  but  the  last  layer.  Thus  the  response 
calculation  may  alternate  between  either  choice  of  independent 
Airy  functions  when  convenient.  If  the  travelling  wave  choice, 

Ai  (-ze_  “^7T//^)  ,  had  been  used,  the  propagator  in  algorithms 
II  and  III  would  be  dominated  by  exponentially  large  terms 
of  nearly  equal  magnitude  and  opposite  sign  whenever  the  para¬ 
meter  N  defined  in  the  Appendix  A  equaled  si.  In  this  case  the 
travelling-standing  choice  avoids  the  numerical  difficulties 
of  the  travelling  choice. 

Algorithm  I  always  involves  fewer  mathematical  operations 
than  II  and  III  if  all  of  the  terms  in  II  and  III  are  cal¬ 
culated.  Branching  between  the  algorithms  can  be  allowed  at 
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non-zero  values  of  m  ^ ,  maximizing  the  use  of  algorithm  I. 

Care  must  be  taken,  however,  that  the  branching  value  of  m^  is 
small  enough  that  precision  is  not  lost  in  algorithm  I. 

The  exponential  scale  factor  accumulated  from  each  redefi¬ 
nition  of  Y^  cancels  in  the  numerator  anc  denominator  of  the 
integrand  until  the  boundary  below  the  source  layer  is 
reached.  The  scale  factor  for  propagation  below  the  source  bound¬ 
ary  can  thus  be  discarded.  The  scale  factors  separately  ac¬ 
cumulated  by  the  numerator  and  denominator  for  propagation  through 
and  above  the  source  layer  can  be  saved  and  used  to  scale  the 
final  result  for  the  integrand. 

Only  tunneled  energy  can  reside  in  a  layer  at  ray  para¬ 
meters  corresponding  to  rays  bottoming  far  above  the  layer. 

These  p  points  occur  whenever  the  parameter  N  in  the  Appendix 
equals  ±1.  The  parameters  |  u>t  j  and  j  u-t  c  |  defined  in  Ap¬ 
pendix  A  then  serve  to  quantify  the  tunneling  phenomenon.  When 
the  tunneled  energy  in  a  layer  is  sufficiently  small  for  both 
P  and  SV  waves  (  |wt^>10  and  }  lut  ^  |  >10  ),  the  propagator  for 
the  layer  can  be  ignored  and  the  matrix  redefinitions  begun 
in  the  next  higher  layer. 


EXAMPLE  USES  OF  GENERAL  SOLUTION  FORM 


I.  Inhomogeneous  Unlayered  Sphere 


As  an  example  of  the  generality  of  (31) ,  consider  an 
explosive  point  source  at  radius  rg  in  an  unlevered  inhomogeneous 
sphere  with  radius  r0.  For  the  radial  displacement  component 
calculate 


-  j  m 

I  U (r_) P  •  r  =  T  - 

m=-n  n  m=-n  RTXXXz 


(41) 


where  for  an  unlayered  sphere 


m 


l  GjP" *r  =  [l,0,0,0]sr  and  R~!  =  FTt(rc). 


-  j 


m=-n 


3-3 


(42) 
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Next  evaluate  the  source  vector  sr  from  the  radiation 
factors  defined  in  Appendix  B.  For  an  explosive  point  source 
the  moment  tensor  is  isotropic ,  giving  in  equation  (BIO) 


In  the  far  field,  the  non-zero  radiation  factors  then  become, 
from  equation  (B9) 


oT~  VcosA) 

s 


Substituting  using  equations  (39)  and  F^j(r0)  using 
the  definition  of  the  inverse  fundamental  matrix  in  equation  (10b) 


gives 


u  (Ar«)  =  / 


-i«p(£)2[^r  '  ^9(3)(rs,Pn(cosA) 
^Vo^K^r  *  g(3)  (r0) 

S  ^  ^  0  P 


dp  (45) 


for  the  radial  component  of  the  Fourier  time  transform  of  dis¬ 
placement.  (n  =  u>p  -1/2) £,  and  n  are  evaluated  at  r„.  The 
denominator  factor, 


can  be  recognized  as  a  form  of  Rayleigh's  dispersion  function 


vith  generalized  vertical  slownesses  substituting  for  the  more 


familiar  plane  wave  cosines. 

This  solution  must  possess  a  ray-mode  duality  equivalent 
to  that  of  the  layered  half  space  described  by  Pekeris  (1948) . 
To  demonstrate  this  duality,  reduce  (45)  to  either  (1)  a 
sum  of  modes  of  free  oscillation  of  the  sphere  or  (2)  to  a  sum 
of  individual  displacements  due  to  all  possible  seismic  rays 


interacting  with  the  free  surface. 
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A  modal  representation  can  be  derived  by  omitting  the 

Poisson  transform  of  the  n  summation  and  evaluating  the  inverse 

Fourier  transform  by  contour  integration,  giving  a  double  summation 

over  radial  order  number  n  and  poles  .  cc  of  Rayleich's  function 

u  n 

in  complex  frequency  (e.g.,  Sato  et  al.,  1963).  An  alternate 
mode  representation  can  be  derived  by  the  spectral  method,  which 
evaluates  the  p  integral  in  (45)  at  the  discrete  frequency  points 
needed  for  the  application  of  a  fast  Fourier  transform.  The 
p  integral  is  then  evaluated  by  contour  integration,  which 
together  with  the  inverse  FFT  results  in  a  double  summation 
over  frequency  a>n  and  poles  ^pn  in  the  complex  p  plane. 

The  validity  of  ray  theory  improves  with  increasing 
frequency.  A  ray  representation  can  thus  be  most  simply 
illustrated  by  incorporating  a  high  frequency  approximation  of 
the  radial  eigenfunctions  in  (45) .  At  sufficiently  high 
frequency  and  small  ray  parameter,  the  first  term  in  an 
asymptotic  representation  of  the  Airy  functions  can  be  used  to 
approximate  the  radial  eigenfunctions.  This  approximation, 
which  is  equivalent  to  the  WKBJ  approximation,  can  be  used  to 
expand  the  integrand  of  (45)  into  a  sum  of  phasor  terms  that 
can  be  identified  with  seismic  rays.  The  displacement  gen¬ 
erated  by  each  ray  can  then  be  determined  by  either  the 
spectral  method  combined  with  saddle  point  integration  (Richards, 
1973),  or  by  Chapman's  (1978a)  technique  of  WKBJ  seismograms. 

II.  Inhomogeneous  Layered  Sphere 

Again  consider  a  buried  explosive  point  source  but  now 
evaluate  the  matrix  products  of  equation  (31)  .  As  a  specific 
example  take  the  layered  earth  model  shown  in  Figure  2,  and 
an  explosion  source  at  1  km.  depth.  Radial  displacement  is 
calculated  by  the  spectral  method  using  a  real  fast  Fourier 
transform: 

ur(t  ,♦,*>=  Re /Y^/ 

r  0  f  Dp_  CT7 


dpdu) 


(46) 
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where  F  is  a  contour  in  the  complex  p  plane. 

The  matrix  products  in  (46)  are  evaluated  at  points  in 
the  complex  ray  parameter  plane,  applying  the  methods  of 
Abo-Zena  (1979)  as  previously  described  in  the  section 
EVALUATION  OF  PROPAGATOR  PRODUCTS .  The  approximations  to 
the  radial  eigenfunctions  at  each  layer  are  constructed  to  be 
continuous  functions  in  the  complex  ray  parameter  plane 
(see  Appendix  A),  allowing  the  integration  path  to  be  deformed 
off  the  real  p  axis.  Deformation  of  the  p  contour  below  the 
real  p  axis  reduces  the  large  oscillations  of  the  integrand 
associated  with  poles  on  and  near  the  real  p  axis.  Numerical 
integration  thus  becomes  a  computationally  efficient  procedure 
without  invoking  an  unreasonably  attenuating  earth  model  to 
move  the  p  poles  off  the  real  axis. 

Residues  at  poles  in  the  first  quadrant  of  the  p  plane 
(Figure  3)  describe  seismic  radiation  at  the  free  surface. 
Application  of  the  reflectivity  boundary  condition  to  the  lower¬ 
most  layer  boundary  produces  poles  in  the  fourth  quadrant  whose 
residues  describe  transmitted  radiation  into  the  last  layer. 

The  integration  contour  in  the  p  plane  must  therefore  be 
deformed  such  that  these  fourth  quadrant  poles  are  excluded 
from  the  response  at  the  free  surface.  Application  of  the 
decay  boundary  condition  in  the  lowermost  layer  eliminates 
the  fourth  quadrant  poles  at  the  expense  of  packing  poles  more 
densely  along  the  real  ray  parameter  axis  at  small  ray  para- 
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meter  values.  The  reflectivity  boundary  condition  is  thus  pre¬ 
ferred  for  the  deformation  of  the  left  end  of  the  contour 
across  the  real  ray  parameter  axis. 

Provided  that  the  Legendre  function  P^fcosl)  is  not  decom¬ 
posed  into  horizontal  travelling  waves  as  in  equation  (Bli) , 
the  integrand  vanishes  as  p  goes  to  zero  and  the  p  contour  may 
end  at  the  origin  as  in  the  antipode  problem  described  by 
Rial  and  Cormier  (1979) .  In  the  distance  range  for  which  sample 
seismograms  were  synthesized  (75-150  km.), near  vertically 
incident  body  waves,  however,  were  not  important  in  the  complete 
response.  Rays  having  take-off  angles  more  vertical  and  bot¬ 
toming  at  depths  deeper  than  the  direct  P  wave  reach  the 
receiver  only  by  suffering  partial  reflections  from  layer 
boundaries.  Such  arrivals  are  consequently  much  smaller 
in  amplitude  and  delayed  in  time  relative  to  the  earlier  portion 
of  the  waveform.  These  arrivals  can  be  excluded  by  deforming 
the  left  end  of  the  contour  across  the  real  p  axis  into  the 
left  portion  of  the  first  quadrant  where  the  integrand  decays. 
The  p  point  at  .which  the  contour  may  be  so  deformed  can  be 
estimated  from  the  lowest  bottoming  depth  of  the  direct  P 
waves  in  the  distance  range  of  calculation.  The  range  of 
distances,  frequencies,  and  ray  parameters  needed  for  the  syn¬ 
thetics  shown  in  Figures  4-5  required  the  inclusion  of  only 
the  Q  ^  ^  travelling  wave  in  the  decomposition  of  the  Legendre 
function,  and  the  left  end  of  the  integration  contour  was 
deformed  across  the  real  p  axis. 

Decay  of  the  integrand  of  equation  (46)  also  occurs  along 
the  real  p  axis  at  large  values  of  p  (e.g.  ,  Frazer,  1977). 

The  integrand  decays  more  rapidly,  however,  by  deforming  the 
right  end  of  the  contour  upward  into  the  upper  right  portion 
of  the  first  quadrant  around  the  last  ray  parameter  pole 
(Figure  3) .  The  residue  at  this  pole  accounts  for  the  funda¬ 
mental  Rayleigh  wave.  For  the  earth  model  chosen  here  the 
position  of  the  fundamental  Rayleigh  pole  is  well  separated 
in  the  complex  p  plane  from  all  other  poles.  It  is  thus 
convenient  to  evaluate  the  displacement  response  in  two 
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stages  (Figure  4),  separately  calculating  the  response  due  to 
the  fundamental  Rayleigh  pole  with  one  contour  and  the  response 
of  all  other  modes  with  another  contour. 

Figures  4-5  show  unfiltered  traces  of  synthetic  displacements . 
The  contour  deformations  shown  in  Figure  3  do  nor  introduce 
truncation  pulses  common  to  methods  that  apply  a  phase  velocity 
filter  to  the  contour  integration.  The  earliest  disturbance 
in  the  wave-forms  represents  the  direct  P  wave.  Some  numerical 
noise  occurs  at  the  Nyguist  frequency  (2  Hz.). 

The  real  axis  contour  can  be  alternately  deformed  to 
completely  enclose  p  poles  in  the  first  quadrant.  The  steps 
in  the  contour  deformation  would  be  identical  to  those  described 
by  Ugincius  and  Oberall  (1968)  for  the  problem  of  the  elastic 
cylinder.  The  spectral  response  could  then  be  evaluated  by 
applying  Cauchy's  residue  theorem  to  the  p  poles.  The  mode 
sum  obtained  would  be  analogous  to  that  obtained  by  Harvey  (1979) 
for  planar  homogeneous  layers. 

In  either  the  numerical  integration  method  or  the  mode 
sum  method  most  of  the  layer  calculations  can  be  saved  and  catalogued 
for  use  at  different  distances  with  different  sources  at  different 
depths.  To  achieve  this  economy  in  the  numerical  integration 
method,  catalogue  the  d?  element  of  Y  matrix  and  all  elements 
of  the  Ys_^  matrix  at  the  frequency  and  ray  parameter  points 
needed  for  the  FFT  and  p  integration. 

MODEL  SPECIFICATION 

A  spherical  earth  model  need  not  be  first  flattened  before 
the  inhomogeneous  layers  are  specified.  The  zero  order  asymptotic 
solutions  account  for  the  curvature  of  layer  boundaries  as  well 
as  inhomogeneity  in  the  layers.  Formulae  for  a  flattened  model 
are  redundant.  The  expressions  for  t  and  Q  needed  for  the  calcula¬ 
tion  of  Airy  functions  in  Appendix  A  are  ecuivalent  to  those  for 

m 

a  spherical  model;  P  (cosA)  is  approximately  replaced  by 
m  n 

(kx)  J  (kx)  with  the  use  of  Szegos  (1934)  asvmptotic  expression 
m 

for  Pm.  (See  for  example  Chapman,  1973,  or  Muller,  1977.) 

n 

The  asymptotic  solution  used  in  this  paoer  restricts  the 
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velocity  profiles  to  having  only  one  turning  point.  The  presence 
of  a  low  velocity  zone  would  thus  preclude  an  earth  model  having 
only  one  layer  (a  simple  inhomogeneous  sphere) . 

How  accurately  do  the  zero  order  asymptotics  model  the 
radial  eigenfunctions  of  an  inhomogeneous  layer?  The  accuracy 
depends  on  both  frequency  and  the  magnitude  of  boundary  curva¬ 
ture  and  first  and  higher  order  radial  derivatives  of  density 
and  elastic  moduli.  Higher  frequency  allows  larger  magnitudes 
for  curvature  and  radial  derivatives.  For  a  given  velocity 
profile  the  error  can  be  estimated  from  either  the  higher  order 
terms  in  the  radial  component  of  the  potential  equations 
(Richards,  1976) ,  or  from  the  first  higher  order  solution  term  for 
the  fundamental  matrix  obtained  by  Woodhouse  (1978) .  with  the 
use  of  smooth  relatively  linear  velocity  gradients,  the  zeroth 
order  asymptotics  can  be  used  in  the  frequency  band  (0.2  -  10  Hz.) 
with  a  relative  error  in  a  radial  eigenfunction  bounded  by  0.1% 
at  0.2  Hz.  for  a  velocity  gradient  of  0.1  sec  . 

Asymptotic  series  type  solutions,  however,  fail  to  properly 
account  for  narrow  angle  reflections  from  regions  having  rapid 
velocity  variations  over  short  distances  (Chapman,  1978b) . 

Layer  boundaries  should  thus  be  introduced  at  points  best 
described  by  a  first  order  discontinuity  such  as  the  moho.  If 
the  errors  in  the  zeroth  order  solution  become  unacceptably 
large  in  a  given  frequency  band,  then  the  layer  should  be  broken 
up  into  thinner  layers  of  weaker  velocity  and  density  gradient. 
Inflections  and  kinks  in  the  velocity  profile,  expressed  by 
a  large  value  of  second  and  higher  order  radial  derivatives, 
should  be  avoided.  In  many  practical  applications  at  500  km. 
or  less  distance,  no  more  than  four  to  five  inhomogeneous  layers 
should  be  necessary  to  describe  the  velocity  variation  from 
crustal  layers,  high  velocity  mantle,  and  low  velocity  zone. 

CONCLUSIONS 


The  preceding  has  demonstrated  the  practicality  of  including 
an  inhomogeneously  layered  model  in  a  calculation  of  the  complete 
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displacement  response  at  short  distances.  The  isolation  of  the 
fundamental  surface  wave  pole  in  the  complex  p  plane  allows  the 
synthesis,  if  desired,  of  only  the  fundamental  surface  wave  mode 
by  numerical  integration. 

When  only  a  small  subset  of  body  waves  at  large  distances 
is  desired,  the  inhomogeneous  layer  matrices  and  response  function 
calculation  can  be  modified  for  a  reflectivity  type  calculation 
analogous  to  that  of  Fuchs  and  Muller  (1971)  for  planar  homogeneous  lay¬ 
ers.  For  such  subsets  of  body  waves  it  would  not  be  difficult  to 
extend  the  decomposition  and  multiplication  scheme  used  for  the 
evaluation  of  the  layer  matrix  to  the  boundary  matrix 
used  in  the  reflection  matrix  method  described  by  Kennet  (1974) , 
where 


Si  =  H-I<ri-I>£i<ri-1> 

In  this  method  the  fundamental  matrices  F.^  and  must  be 
defined  such  that  they  can  be  partitioned  into  four  2X2  matrices 
each  of  which  can  be  associated  solely  with  up-  or  down-going 
waves. 

In  any  of  these  applications  the  layer  inhomogeneity 
minimizes  the  number  of  layers  needed  to  describe  the  response 
function  and  thereby  the  computation  time  required  to  synthesize 
displacement.  The  zero  order  asymptotics  adequately  model  the 
behavior  of  the  radial  eigenfunctions  of  an  inhomogeneous  layer 
in  the  frequency  band  (0.2  -  10  Hz.)  for  the  velocity  gradients 
common  to  wide  regions  of  the  crust  and  upper  mantle. 
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Figure  1 
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Figure  4 


Figure  5 


FIGURE  CAPTIONS 


Layering  index  scheme. 

Earth  model  based  on  the  southern  California  crustal 
model  of  Kanamori  and  Hadley  (1975)  and  model  CIT  109 
of  the  upper  mantle  (Archambeau  et  al.,  1969). 

Paths  in  the  complex  ray  parameter  plane  for  the 
evaluation  of  displacement  in  a  spherical  earth  model 
having  a  single  discontinuity  at  which  the  P  and  S 
velocity  increase  discontinuously .  Position  of  poles 
in  the  P-SV  response  are  shown  schematically  by  x's. 
Rayleigh  poles  lie  near  the  real  axis.  Franz  poles 
emanate  upward  and  downward  from  the  real  p  axis  at 
approximately  ±60  degree  angles.  The  poles  of  the 
P-SV  response  of  the  earth  model  shown  in  Figure  2 
differ  from  this  form  only  in  having  additional 
Franz  poles  in  the  first  quadrant  associated  with 
higher  layer  boundaries. 

Displacement  calculated  from  the  spectral  response 
(0  -  2  Hz.)  given  by  the  numerical  integration  along 
the  contours  shown  in  Figure  3.  Polarity  is  reversed. 
Source  was  an  explosion  at  1  km.  depth.  Time  domain  re¬ 
cords  were  obtained  by  inverse  FFT  using  128  frequencies. 
Synthetic  displacement  at  75  to  150  km. 
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APPENDIX  A:  RADIAL  EIGENFUNCTIONS 

Practical  Evaluation  of  Airy  Functions 

The  radial  eigenfunctions  are  defined  in  terr_s  of  the 
Airy  function  Ai : 
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The  exact  radial  eigenfunctions  of  an  inhomrgeneous  layer 
are  analytic  functions  in  the  complex  p  plane.  The  fractional 
powers  used  in  defining  the  arguments  cf  the  Airy  functions, 
however,  place  branch  cuts  in  the  complex  p  plane  cf  the  asymptotic 
approximations  of  the  radial  eigenfunctions.  These  branch  cuts 


may  be  removed  by  multiplying  wt  by  e‘'27TN  ,  where  N=0,  ±1 
depending  on  the  phase' of  the  complex  t  and  ray  parameters  p. 
Cormier  (1976a)  discusses  how  N  may  be  determined  in  a  complex 
velocity  profile  when  subroutine  packages  for  Hankel  functions 
or  order  1/3  calculate  the  Airy  functions.  The  Airy  functions, 
however,  can  be  most  efficiently  calculated  with  the  exact  and 
asymptotic  formulae  given  in  Abramowitz  and  Stegun  (1964) . 

P.G.  Richards  (personal  communication)  has  used  this  approach, 
patching  the  exact  and  approximate  formulae  together  when 
|  d)T  |  =  3.5.  The  subroutine  LANGER  written  and  tested  by  Richards 
together  with  an  improved  algorithm  for  calculating  N  (a  portion 
of  his  subroutine  ROMTAU)  have  been  used  in  the  calculations  of 
this  paper.  The  integers  rru  needed  in  propagator  multiplication 
scheme  of  equations  (39-40)  can  be  constructed  from  normal  LANGER 
output. 


Practical  Evaluation  of  the  Delay  Time  t 

The  integral  defining  t  may  be  (1)  analytically  evaluated 
using  the  velocity  profile  v  =  ar  in  each  layer  (Richards,  1973, 
1976;  Choy,  1977;  Frazer,  1977;  Mondt,  1977);  (2)  numerically 

integrated  in  an  analytic  velocity  profile  given  by  a  polynomial 
in  radius  (Cormier  and  Richards,  1977);  or  (3)  analytically 
integrated  by  parameterizing  the  logarithm  of  radius  in  the 
Bullen  parameter  (In  =  r/v(r))  (Woodhouse,  1974). 

Method  (3)  is  equivalent  to  the  parameterization  described 
by  Cerveny  et  al.  (1977)  and  used  by  Germany,  et  al.  (1979)  for 
flat  inhomogeneous  layers.  It  is  usually  the  most  computationally 
efficient.  For  the  calculation  of  the  seismograms  in  Figures  4-5, 
each  layer  of  the  velocity  model  in  Figure  2  was  parameterized 
using  method  (3)  with 

In  (r/r  0 )  *  a  +  b(n/ri0)2  (A2) 

where  r0  and  n0  are  the  values  of  r  and  n  at  the  top  of  the  layer. 
The  constants  a  and  b  were  determined  from  fhe  values  of  radius 
and  velocity  at  the  boundaries  of  the  layer.  The  delay  time 


was  then  given  by 


T(p)  =  2n0  [1  -  (p/rj0 )  2  ]  3/2  b/3 


(A3) 


Incorporation  of  Attenuation 

Incorporation  of  attenuation  in  any  of  the  above  methods  for 
computing  t  can  be  accomplished  by  a  complex  velocity  profile 
(e.g.,  Cormier  and  Richards,  1976).  When  a  complex  velocity 
profile  can  be  specified  by  the  form 

O(r)  =  v0 (r) *  <1  +  t)  (A4 ) 

method  (3)  reduces  to  the  form 

tA(p)  =  (1  +  e)_1Tu(p*(l  +  6))  (A5) 

where  e  is  a  small  complex  number  constant  in  radius,  vo (r)  a 
real  velocity,  ta  the  delay  time  in  the  attenuating  layer,  and 
Ty  the  delay  time  in  the  layer  having  the  real  velocity  profile 
evaluated  at  the  ray  parameter  p=p* (1+6) .  To  include  the  dispersive 
effect  of  linear  attenuation  mechanisms  £  must  have  non-zero 
real  and  imaginary  parts,  each  of  which  is  frequency  dependent 
(Futterman,  1962) .  If  dispersion  can  be  ignored  in  the  frequency 
band  of  interest,  then  £  may  be  taken  to  be  pure  imaginary  and 
frequency  independent. 
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emphasize  the  major  properties  of  the  source  radiation  pattern  and 
are  defined  by 
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Specification  of  a  Seismic  Point  Source 


The  expressions  in  equations  (B2)  can  be  simplified  by 
specifying  the  seismic  source  vectors  y  and  Y ' ♦  Consider  a 
seismic  source  specified  by  the  Fourier  time  transform  f 
of  a  vector  £  of  body  force  density  acting  at  a  point.  As  in 
Saito  (1967)  ,  y  and  y1  are  determined  from  an  expansion  of 
f  in  vector  spherical  harmonics: 
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* denoting  the  compex  conjugate. 

f  can  then  be  described  with  the  zeroth  order  (in  frequency) 
moment  rate  tensor  M  by 


f  =  -M* V6 (r  -  r£) 


(B5) 
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Gilbert  (1971),  where  6(r  -  r  )  is  a  three  dimensional  delta 

—  — s 

function  given  at  a  point  with  radius  vector  r  .  In  spherical 
polar  coordinates, 

6(r  -  r  )  =  6(r-rs)6(e-60)6(<D-<^D) 

—  — s  - 

r2sinS 


A  distributed  source  can  be  represented  by  a  point  source 
with  the  device  of  multipolar  expansions  (Archambeau,  1968) . 

A  sufficiently  accurate  approximation  of  the  radiation  pattern 
of  a  distributed  source  a  high  frequency,  however,  can  require 
the  inclusion  of  many  terms  in  such  an  expansion.  The  equivalent 
definition  with  a  moment  tensor  requires  many  higher  order  moment 
tensors  (Gilbert,  1971)  together  with  higher  order  spatial 
derivatives  of  the  delta  function  (Hudson,  1969)  .  The  radiation 
factors  given  here  consider  only  the  first  term  ir.  such  an 
expansion,  a  source  representation  always  valid  at  sufficiently 
low  frequency.  The  same  solution  procedure  may  be  used  when 
higher  order  moment  tensers  are  included. 

The  properties  of  M  have  been  discussed  by  Gilbert  (1971) . 

By  conservation  of  angular  momentum,  M  must  be  a  symmetric 
matrix  (Mre  =  Mgr,  =  K^)  .  For  a  point 

double  couple,  quadrapole  source  M  must  have  zero  trace 
(Mrr  +  =  0)  .  For  a  pure  explosive  point  source, 

M  must  be  isotropic,  i.e.,  all  off  diagonal  elements  are  equal 
(Mr6  =  Mr^  =  Mg^  *  0) ,  and  all  diagonal  elements  are  equal 
(Mrr  *  Mgg  =  =  M0).  An  exposion  having  some  tectonic  stress 

release  may  be  represented  by  an  M  in  which  the  off  diagonal 
elements  are  non-zero.  The  elements  of  M  have  been  defined  by 
Gilbert  (1971)  and  Aki  and  Richards  (1980)  in  terms  of  earth¬ 
quake  fault  plane  parameters. 

Results  for  the  Source  Radiation  Vectors  sr,  etc. 

Simplified  expressions  for  js  etc.  may  now  be  found  from 


substituting  equations  (10b) , (B3) , (B4) , (B5) ,  ana  (B6)  into 
equations  (B2)  and  applying  the  addition  theorem  for  vector 
spherical  harmonics.  For  examples  of  the  simplification 
procedure  see  Singh  and  Ben-Menahem  (1969)  and  Ward  (1979) . 
In  the  results  that  follow,  L  is  the  epicentral  distance. 

The  azimuthal  angles  C  and  E  are  defined  by  Singh  and  Ben- 
Menahem  for  a  generalized  co-ordinate  system  at  the  source 
and  by  Richards  (1978)  and  Ward  (1979)  for  a  north-east-down 
(NED)  co-ordinate  system. 


c(3),  c(1),  h^,  and  h^  are  evaluated  at  the  source  radius 


Ill 


rc  in  the  analytic  velocity  profiles  of  the  source  layer. 

®  P 

The  radiation  factors  consist  of  far  field  terms  ( 
SV'  £  SYL1  anc^  near  field  terms  (  Xpr  X-SV'  : 


(B8) 


with 

FIp  =  [p2/r2LF  +  l/a2LF]Pn (cosA)  -  2£p/r  (cosA) /( iwp) 

FJp  =  [p2/r2  LF  +  l/a|  L^]  P^  (cosA)  +  2£p/r  L*P^(cosA) /Uup) 
F£Sv  ~  (Qg  “  P2/r2)  L^Vp^(cosA)  /(iwp)  -  np/r  L^P^cosA) 

^£sv  =  (Q2  -  p2/r2)L^Vp^(cosA)/(£u!p)  +  fip/r  L^VPn(cosA) 

(B9) 

^[>SH  =  nL^HPn(cosA)  -  p/r  Lg^P^  (cosA)  / (iup) 

FPcn  *  -nL1  P  (cosA)  -  p/r  L2  P*  (cosA)  /(iup) 

>tSH  SH  n  SH  n 

%  ■  NXp  =  p2/r2  L*  cotA  P^  (cosA) / (u2p2 ) 
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* 


Ni 

NI 


sv 

SH 


Ni 

^SH 


-np/r  L1*  cotAP^ (cosA) /(u2p2 ) 
P  s 

np/rLpgCOtAP^ (cosA) /{w2p2 ) 


(39) 

(cor.t ’d) 


-2nLgHcotAP^ (cosA) /(u2p2)  -  p/r  L2H [2cotA?n (cosA) / (  iup)  + 
(2csc2A  -  4cot2A)  P^  (cosA)/(i  u3p  5)  ] 


2nLggCOtAP^  (cos A)  /  (aj2p  2  )  -  p/r  l|w  [2cotAPn  (cosA)  /  (  icp)  + 
(2csc2A  -  4cot2 A) P^ (cosA) / (iw3p 3) 


where  p  independent  factors  are  given  by : 


= 

Me9c°S23 

+  Me<^sin2 3  +  M^sin 23  -  Mrj. 

LP 

= 

Mr0cos=  • 

+  M  ,sin3 
r4> 

Li 

= 

M 

rr 

Lsv 

as 

LP 

Lsv 

as 

LP 

T  1 

"SH 

as 

,cos3 

r? 

-  Mresin= 

lsk 

= 

sinS (M  . 

<P<? 

-  Keel/2  +  cos2»M6t 

as 

cos23MQe 

+  2M„ , sin23  +  2M, ,sin23 

u  <p  $$ 

(BIO) 


The  radiation  factors  defined  remain  valid  as  A  vanishes.  For 

A  not  near  0,  it  is  possible  to  divide  P  (cosA)  into  travelling 

(1)  (2)  n 
wave  functions  Qn  (cosA)  and  Qn  (cosA) , 

Pn(COsA)  =  Qn^(cosA)  +  Qn^(cosA)  (Bll) 

For  large  u)p  these  behave  asymptotically  as 
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(1)  _ 

,  <2>  -  ./  1.  ±i[u)P-i:/4] 

n  \  2wpTTsinA 


(B12 ) 


l  SP 

Nussenzweig  (1965)  .  With  P  =  and  ecuations  (B11-B12)  , 

n  cl. 

it  follows  that 

P^(cosA)  =  iojpPn  (cosA)  .  (B13) 


Thus,  for  large  topA  the  far  field  terms  behave  like  P  (cosA) ,  and 

n  2 

the  near  field  terms  like  Pn (cosA) / (wpA)  and  Pn (cosA) / (wpA)  . 

The  radiation  factors  in  equations  (BIO)  are  essentially 
the  same  as  those  given  in  Aki  and  Richards  (1980)  and  Ward  (1979) . 
In  the  far  field  only  three  p  independent  factors  are  necessary 
to  describe  the  radiation  of  P  waves,  and  two  p  independent 
factors  to  describe  the  radiation  of  SV  or  SK  waves.  In  order 
to  fully  account  for  the  radiation  pattern  of  the  source  in  the 
near  field,  however,  the  l£s  factor  must  be  added  to  the  results 
reported  in  Aki  and  Richards  and  the  Legendre  functions  not 
divided  into  travelling  wave  functions  near  A=0. 


